{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.9" }, "colab": { "name": "403_Adams Moulton Example.ipynb", "provenance": [], "include_colab_link": true } }, "cells": [ { "cell_type": "markdown", "metadata": { "id": "view-in-github", "colab_type": "text" }, "source": [ "\"Open" ] }, { "cell_type": "markdown", "metadata": { "id": "oqtXx_ihb5Sp" }, "source": [ "# Adams Moulton\n", "#### John S Butler \n", "john.s.butler@tudublin.ie \n", "[Course Notes](https://johnsbutler.netlify.com/files/Teaching/Numerical_Analysis_for_Differential_Equations.pdf) [Github](https://github.com/john-s-butler-dit/Numerical-Analysis-Python)\n", "\n", "\n", "The Adams Moulton method is an implicit multistep method. This notebook illustrates the 2 step Adams Moulton method for a linear initial value problem of the form\n", "\\begin{equation} y^{'}=t-y, \\ \\ (0 \\leq t \\leq 2)\\end{equation}\n", "with the initial condition\n", "\\begin{equation}y(0)=1.\\end{equation}\n", "The video below walks through the notebook." ] }, { "cell_type": "code", "metadata": { "id": "MDJbhHZ6b5Sr", "colab": { "base_uri": "https://localhost:8080/", "height": 336 }, "outputId": "bddd23e8-f464-4b63-958d-365ef387dabc" }, "source": [ "from IPython.display import HTML\n", "HTML('')\n" ], "execution_count": 1, "outputs": [ { "output_type": "execute_result", "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "execution_count": 1 } ] }, { "cell_type": "markdown", "metadata": { "id": "YwzVuSStb5Sw" }, "source": [ "## Python Libraries" ] }, { "cell_type": "code", "metadata": { "id": "URyGrrhKb5Sx" }, "source": [ "import numpy as np\n", "import math \n", "import pandas as pd\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt # side-stepping mpl backend\n", "import matplotlib.gridspec as gridspec # subplots\n", "import warnings\n", "\n" ], "execution_count": 2, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "dYFIAsGRb5S0" }, "source": [ "### Defining the function\n", "\\begin{equation} f(t,y)=t-y.\\end{equation}" ] }, { "cell_type": "code", "metadata": { "id": "2uv4HlyXb5S1" }, "source": [ "def myfun_ty(t,y):\n", " return t-y" ], "execution_count": 3, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "PO_T_dn9b5S4" }, "source": [ "## Discrete Interval\n", "Defining the step size $h$ from the interval range $a\\leq t \\leq b$ and number of steps $N$ \n", "\\begin{equation}h=\\frac{b−a}{N}.\\end{equation}\n", " \n", "This gives the discrete time steps,\n", "\\begin{equation}t_i=t_0+ih,\\end{equation}\n", "where $t_0=a.$\n", "\n", "Here the interval is $0≤t≤2$ and number of step 4 \n", "\\begin{equation}h=\\frac{2−0}{4}=0.5.\\end{equation}\n", " \n", "This gives the discrete time steps,\n", "\\begin{equation}t_i=0+i0.5,\\end{equation}\n", "for $i=0,1,⋯,4.$" ] }, { "cell_type": "code", "metadata": { "id": "ppaw5pK_b5S4", "colab": { "base_uri": "https://localhost:8080/", "height": 298 }, "outputId": "17afa13a-3ca0-42cb-83bd-5d9c0baf8da9" }, "source": [ "# Start and end of interval\n", "b=2\n", "a=0\n", "# Step size\n", "N=4\n", "h=(b-a)/(N)\n", "t=np.arange(a,b+h,h)\n", "fig = plt.figure(figsize=(10,4))\n", "plt.plot(t,0*t,'o:',color='red')\n", "plt.xlim((0,2))\n", "plt.title('Illustration of discrete time points for h=%s'%(h))" ], "execution_count": 4, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "Text(0.5, 1.0, 'Illustration of discrete time points for h=0.5')" ] }, "metadata": {}, "execution_count": 4 }, { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmwAAAEICAYAAADiGKj0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAedklEQVR4nO3de7xldV3/8ddbRjHEuItcZ1DoAmVeTqBlNipysWQw6fcbQ0XUJgos62eFUl5ACs0C/akZKklAgJHmUCEXEfr9QIkzRAoqMozogCCXGS4jCqKf/ljrOHvO7HPOPu4zc9Y5vJ6Px3nM+q71XWt913d/997vsy5zUlVIkiSpux432w2QJEnS5AxskiRJHWdgkyRJ6jgDmyRJUscZ2CRJkjrOwCZJktRxBjapleS1Sf5/T7mS7D2bbZpIkg8n+fNZ2O/vJvl2knVJdhig/q1JDmyn35rko5u+lZtOkiOTXDLb7ZhMkj3b12eLzbCvlydZ3e7vWTOwvSuSvGEm2ibNNwY2Pab1BopNtP2PJ3nXkNvYIEgCVNUxVXXScK2bdjseD/wNcFBVbV1V905n/ar6i6rabF/GSRa1oXvBTK1fVedU1UEz18qZV1XfbF+fH0xVd9g+At4LHNfu779+zG3MuCTPTLIiyUPtv8+cpO4VSb7Xhs51SW7anG2VBmVgk2bREF+Us2Fn4InAjbPdEJhzfTdfLeTHHA+b6gxgkicAnwbOBrYDzgQ+3c6fyFjo3LqqfnpTtEsaloFNGsD4SzW9Z73SODXJXUkeSPKlJD+XZBlwJPAn7W/uF7b1b03yp0m+CHwnyYIkxye5JcmDSb6c5OVt3Z8FPgw8r93Gfe38Dc7cJfntJCuTrEmyPMmuPcsqyTFJbk5yX5IPJskEx7llktOSfKv9Oa2d91PA2JmH+5JcPsH6r07yjST3Jjlh3LJ3JDm7nX5ikrPbevcluTbJzu2y7ZP8fbv/tUn+pZ2/OMltbd/dCfx9ksf19N29ST6RZPt2l//R0951SZ7Xbud1Sb7SbvviJAsneNk3Wn/82c62b3+v7dsHk5yU5OlJrm7Hwid6g0KSX09yfXvMVyd5xgT7Htv27ydZleSeJH+V5HHtsscl+bO2r+9K8g9JtmmXbXDWrB27JyW5qm3jJUl2nOQY905yZZL72/2e36dtWyZZB2wB/HeSW9r5P9vu774kNyY5rGedjyf52yT/nuQ7wAsnOPSFE7R1UIuBBcBpVfVwVb0fCPCiaW5H6hQDmzS8g4AXAD8FbAP8L+DeqjodOAd4T/ub+8t61nkl8GvAtlX1KHAL8Cvt+u8Ezk6yS1V9BTgG+Hy7jW3H7zzJi4C/bPe7C/AN4Lxx1X4d+EXgGW29gyc4lhOA5wLPBH4B2B/4s6r6GrBfW2fbqtroyy/JvsDfAq8GdgV2AHafYD9Htce6R1vvGOC77bKzgK3a/T0FOLVnvacC29Oc2VkGvBE4HPjVdp9rgQ+2dV/Q096tq+rzSZYAbwV+A9gJ+H/AuRO0caP1J6h3MPAcmn77E+B04FXtsf0czWtNmnu8zgB+pz3mvwOWJ9lygu0CvBwYAZ4NLAFe185/bfvzQuBpwNbABybZzm8BR9P05xOAN09yjCcBl9Ccndod+L/jN9YGoa3b4i9U1dPTXDK/sF33KTSvzTlJes9Y/RZwMvBkYIPL/AO0lTYITvRzfFttP+CLteHfXfwi68dvP3/ZhtOrkiyepJ40awxs0vC+T/MF9DNAquorVXXHFOu8v6pWV9V3Aarqn6rqW1X1w6o6H7iZJiwN4kjgjKq6rqoeBt5Cc0ZuUU+dU6rqvqr6JvA5mkA20bZOrKq7qupumvD46gHbcQTwr1X1H207/hz44QR1v08TWvauqh9U1YqqeiDJLsChwDFVtbaqvl9VV/as90Pg7W1g+C5N0Duhqm5r9/kO4IhMfLn0GOAv29foUeAvgGdOcpZtEO+pqgeq6kbgBuCSqlpVVfcDFwFjN+MvA/6uqq5pj/lM4GGaoDeRd1fVmvZ1O402/NG8Tn/T7mcdzWu+dJLj/vuq+lrbZ59g4tcfmtdmIbBrVX2vqiYKVuM9lyY4nlJVj1TV5cC/9rQZ4NNVdVU7zr833bZW1baT/JzSVtsauH/cNu+neY/286c0oXc3mrB9YZKnD3jM0mZjYJOG1H4xfYDmzM5dSU5P8pNTrLa6t5DkNT2Xyu6jOTMz6KWgXWnOqo21Zx1wL80X0Jg7e6YfovlSm3Jb7fSuE9Ttt+6PjquqvtO2o5+zgIuB89pLn+9pz9DsAaypqrUTrHf3uC/6hcCnevrtK8APaO6362ch8L6e+mtoLpftNkH9QXy7Z/q7fcpjfb0Q+D+9Z4Vojney/u0dJ72vRb/XaQETH/egrz80ZwkD/Gd7WfN1k9TttSuwuqp6Q/o32LBvVzO16bS1n3XA+PffTwIP9qvcBugH218CzgSuAl46zX1Km5yBTRrMd2gu0415au/Cqnp/VT0H2Jfm0ugfjy2aYHs/mt+e3fkIcBywQ3vZ8waaL83JtjHmWzRhYGx7T6I5e3X7FOtNuS1gz3beIO6gCSBj7diqbcdG2jNn76yqfYFforlk+xqaL/Ttk2x06Xds1XHl1cCh4860PLGqbu9Td6z+74yr/xNVdfUA+xrWauDkcfveqqomuiQLPf3Jhq9Fv9fpUTYMi4PY6Bir6s6q+u2q2pXm8u2HMth/b/MtYI+x++x62tU7Dofq06x/krPfz1vbajcCz0g2uE/zGQz+cESx/r0ndYaBTRrM9cBvJNmq/fJ6/diCJL+Y5ID2DNF3gO+x/lLgt2kut0zmSTRfEne32zua5gzbmG8Du2fip9zOBY5O818ZbElzme+aqrp1OgfYs60/S7JTe7P322iethvEBcCvJ3l+29YTmeAzJskLk/x8micFH6C5DPfD9lLyRTQhYbskj0/ygn7baH0YOHnskmbb7iXtsrtpXoenjav/liT7tfW3SfKbE2y73/rD+AhwTDtWkuRJSX4tyUSX6gD+uO2HPYA/AMYeADgX+MMkeyXZmuY1P7+9zDsdGx1jkt9MMnbv4VqasTnRpe1e19CcEfuT9nVbDLyMje+n/LH1PMnZ7+cv2mpX0Jxl/f324Yjj2vkbPSiTZNskB6d5CGZBkiNp7uv7zEy1WZopBjZpMKcCj9CEpzNpHiYY85M0X8ZraS4B3Qv8VbvsY8C+7SWwf+m34ar6MvDXwOfb7f88zWWZMZfTnB24M8k9fda/jOZ+sX+mOcv1dGDpj3WU8C5glOYm7S8B17XzptTew3Us8I9tO9YCt01Q/ak0Ae8BmsuYV9JcJoXmnrnvA18F7gLeNMlu3wcsBy5J8iDwBeCAtj0P0dzgflXb/8+tqk8B76a5FPsAzZnMQyc4no3Wn7ITJlFVo8Bv01w+XwuspHlwYDKfBlbQ/MLwbzTjCZqHF86iecrz6zS/JLzxx2hTv2P8ReCaNE+BLgf+oKpWDbCtR2gC2qHAPcCHgNdU1Ven265htO04nOaM7X00D2oc3s4f+w+cL2qrP55mfN/dtvmNbd2vbc42S4PIhg/SSJK6IEkB+1TVytlui6TZ5xk2SZKkjjOwSZIkdZyXRCVJkjrOM2ySJEkdNyf/ePKOO+5YixYtmu1mSJIkTWnFihX3VNVOw2xjTga2RYsWMTo6OtvNkCRJmlKSb0xda3JeEpUkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeq4GQlsSQ5JclOSlUmO77N8yyTnt8uvSbJo3PI9k6xL8uaZaI8kSdJ8MnRgS7IF8EHgUGBf4JVJ9h1X7fXA2qraGzgVePe45X8DXDRsWyRJkuajmTjDtj+wsqpWVdUjwHnAknF1lgBnttMXAC9OEoAkhwNfB26cgbZIkiTNOzMR2HYDVveUb2vn9a1TVY8C9wM7JNka+FPgnVPtJMmyJKNJRu++++4ZaLYkSdLcMNsPHbwDOLWq1k1VsapOr6qRqhrZaaedNn3LJEmSOmLBDGzjdmCPnvLu7bx+dW5LsgDYBrgXOAA4Isl7gG2BHyb5XlV9YAbaJUmSNC/MRGC7FtgnyV40wWwp8Fvj6iwHjgI+DxwBXF5VBfzKWIUk7wDWGdYkSZI2NHRgq6pHkxwHXAxsAZxRVTcmOREYrarlwMeAs5KsBNbQhDpJkiQNIM2JrrllZGSkRkdHZ7sZkiRJU0qyoqpGhtnGbD90IEmSpCkY2CRJkjrOwCZJktRxBjZJkqSOM7BJkiR1nIFNkiSp4wxskiRJHWdgkyRJ6jgDmyRJUscZ2CRJkjrOwCZJktRxBjZJkqSOM7BJkiR1nIFNkiSp4wxskiRJHWdgkyRJ6jgDmyRJUscZ2CRJkjrOwCZJktRxBjZJkqSOM7BJkiR1nIFNkiSp4wxskiRJHWdgkyRJ6jgDmyRJUscZ2CRJkjrOwCZJktRxBjZJkqSOM7BJkiR1nIFNkiSp42YksCU5JMlNSVYmOb7P8i2TnN8uvybJonb+S5KsSPKl9t8XzUR7JEmS5pOhA1uSLYAPAocC+wKvTLLvuGqvB9ZW1d7AqcC72/n3AC+rqp8HjgLOGrY9kiRJ881MnGHbH1hZVauq6hHgPGDJuDpLgDPb6QuAFydJVf1XVX2rnX8j8BNJtpyBNkmSJM0bMxHYdgNW95Rva+f1rVNVjwL3AzuMq/MK4LqqengG2iRJkjRvLJjtBgAk2Y/mMulBk9RZBiwD2HPPPTdTyyRJkmbfTJxhux3Yo6e8ezuvb50kC4BtgHvb8u7Ap4DXVNUtE+2kqk6vqpGqGtlpp51moNmSJElzw0wEtmuBfZLsleQJwFJg+bg6y2keKgA4Ari8qirJtsC/AcdX1VUz0BZJkqR5Z+jA1t6TdhxwMfAV4BNVdWOSE5Mc1lb7GLBDkpXAHwFj//XHccDewNuSXN/+PGXYNkmSJM0nqarZbsO0jYyM1Ojo6Gw3Q5IkaUpJVlTVyDDb8C8dSJIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcTMS2JIckuSmJCuTHN9n+ZZJzm+XX5NkUc+yt7Tzb0py8EA7XLECFi2Cc86ZieZrPjrnnGaMPO5xjhVNzfGiQTlWNB3teHkOPGfYTS0YdgNJtgA+CLwEuA24NsnyqvpyT7XXA2urau8kS4F3A/87yb7AUmA/YFfgsiQ/VVU/mHLH3/gGLFvWTB955LCHofnknHOasfHQQ03ZsaLJOF40KMeKpmP8eBnSTJxh2x9YWVWrquoR4Dxgybg6S4Az2+kLgBcnSTv/vKp6uKq+DqxstzeYhx6Ct74VFi+Gs89eP2/xYjj//KZ8//1N+ZOfbMr33NOUL7ywKd95Z1P+zGea8urVTfmyy5ryqlVN+corm/JNNzXlq69uyjfc0JSvvbYpX399U77++qZ87bVN+YYbmvLVVzflm25qylde2ZRXrWrKl13WlFevbsqf+UxTvvPOpnzhhU35nnua8ic/2ZTvv78pn39+Ux4bIGef3ZS///2m/PGPN+UxH/kIHHjg+vKHPgSHHrq+/L73wWGHrS+/973wilesL59yCixdur580knwqletL7/tbXD00evLb3nL+g84gDe/GY49dn35TW9qfsYce2xTZ8yyZc02xhx9dLOPMa96VbPO+DfIQw/BCSc0bX/ve9fPP+yw5hjHHHpo0wdjDjyw6aMxixc3fQhNnzr25v7YO+GE/uPlmGOmP/ZOOml9eenSpo1jHHtzf+xNNFaOPXb2P/cce025S2Ov33gZwkwEtt2A1T3l29p5fetU1aPA/cAOA64LQJJlSUaTjG6wYPXqftX1WDb2Rhrvm9/cvO3Q3DDRuFi3bvO2Q9030ViZ6DNHj20z/J2TqhpuA8kRwCFV9Ya2/GrggKo6rqfODW2d29ryLcABwDuAL1TV2e38jwEXVdUFk+1zJKkfpbaFC+HWW4c6Bs0zixY1lyrGc6yoH8eLBuVY0XT0jJcRYLQqw2xuJs6w3Q7s0VPevZ3Xt06SBcA2wL0DrjuxrbaCk0+efos1v518cjM2ejlWNBHHiwblWNF09BsvQ5iJwHYtsE+SvZI8geYhguXj6iwHjmqnjwAur+bU3nJgafsU6V7APsB/DrTXhQvh9NO90VMbO/LIZmwsXAiJY0WTc7xoUI4VTUfveJkBQ18SBUjyUuA0YAvgjKo6OcmJwGhVLU/yROAs4FnAGmBpVa1q1z0BeB3wKPCmqrpoqv2NjIzU6OjoVNUkSZJmXZIVVTUy1DZmIrBtbgY2SZI0V8xEYPMvHUiSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOm6owJZk+ySXJrm5/Xe7Ceod1da5OclR7bytkvxbkq8muTHJKcO0RZIkab4a9gzb8cBnq2of4LNteQNJtgfeDhwA7A+8vSfYvbeqfgZ4FvDLSQ4dsj2SJEnzzrCBbQlwZjt9JnB4nzoHA5dW1ZqqWgtcChxSVQ9V1ecAquoR4Dpg9yHbI0mSNO8MG9h2rqo72uk7gZ371NkNWN1Tvq2d9yNJtgVeRnOWTpIkST0WTFUhyWXAU/ssOqG3UFWVpKbbgCQLgHOB91fVqknqLQOWAey5557T3Y0kSdKcNWVgq6oDJ1qW5NtJdqmqO5LsAtzVp9rtwOKe8u7AFT3l04Gbq+q0KdpxeluXkZGRaQdDSZKkuWrYS6LLgaPa6aOAT/epczFwUJLt2ocNDmrnkeRdwDbAm4ZshyRJ0rw1bGA7BXhJkpuBA9sySUaSfBSgqtYAJwHXtj8nVtWaJLvTXFbdF7guyfVJ3jBkeyRJkuadVM29q4sjIyM1Ojo6282QJEmaUpIVVTUyzDb8SweSJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHXcUIEtyfZJLk1yc/vvdhPUO6qtc3OSo/osX57khmHaIkmSNF8Ne4bteOCzVbUP8Nm2vIEk2wNvBw4A9gfe3hvskvwGsG7IdkiSJM1bwwa2JcCZ7fSZwOF96hwMXFpVa6pqLXApcAhAkq2BPwLeNWQ7JEmS5q1hA9vOVXVHO30nsHOfOrsBq3vKt7XzAE4C/hp4aKodJVmWZDTJ6N133z1EkyVJkuaWBVNVSHIZ8NQ+i07oLVRVJalBd5zkmcDTq+oPkyyaqn5VnQ6cDjAyMjLwfiRJkua6KQNbVR040bIk306yS1XdkWQX4K4+1W4HFveUdweuAJ4HjCS5tW3HU5JcUVWLkSRJ0o8Me0l0OTD21OdRwKf71LkYOCjJdu3DBgcBF1fV31bVrlW1CHg+8DXDmiRJ0saGDWynAC9JcjNwYFsmyUiSjwJU1Rqae9WubX9ObOdJkiRpAKmae7eDjYyM1Ojo6Gw3Q5IkaUpJVlTVyDDb8C8dSJIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjktVzXYbpi3Jg8BNs92OjtkRuGe2G9FB9kt/9kt/9svG7JP+7Jf+7Jf+frqqnjzMBhbMVEs2s5uqamS2G9ElSUbtk43ZL/3ZL/3ZLxuzT/qzX/qzX/pLMjrsNrwkKkmS1HEGNkmSpI6bq4Ht9NluQAfZJ/3ZL/3ZL/3ZLxuzT/qzX/qzX/obul/m5EMHkiRJjyVz9QybJEnSY4aBTZIkqeM6FdiSHJLkpiQrkxzfZ/mWSc5vl1+TZFHPsre0829KcvDmbPemNkC//FGSLyf5YpLPJlnYs+wHSa5vf5Zv3pZvWgP0y2uT3N1z/G/oWXZUkpvbn6M2b8s3nQH65NSe/vhakvt6ls3nsXJGkruS3DDB8iR5f9tvX0zy7J5l83WsTNUnR7Z98aUkVyf5hZ5lt7bzr5+J/66gSwbol8VJ7u95r7ytZ9mk77+5bIB++eOePrmh/TzZvl02L8dLkj2SfK79/r0xyR/0qTNzny1V1YkfYAvgFuBpwBOA/wb2HVfn94APt9NLgfPb6X3b+lsCe7Xb2WK2j2kz9ssLga3a6d8d65e2vG62j2EW++W1wAf6rLs9sKr9d7t2ervZPqbN0Sfj6r8ROGO+j5X22F4APBu4YYLlLwUuAgI8F7hmPo+VAfvkl8aOFTh0rE/a8q3AjrN9DLPUL4uBf+0zf1rvv7n2M1W/jKv7MuDy+T5egF2AZ7fTTwa+1ud7aMY+W7p0hm1/YGVVraqqR4DzgCXj6iwBzmynLwBenCTt/POq6uGq+jqwst3efDBlv1TV56rqobb4BWD3zdzG2TDIeJnIwcClVbWmqtYClwKHbKJ2bk7T7ZNXAudulpbNsqr6D2DNJFWWAP9QjS8A2ybZhfk7Vqbsk6q6uj1meOx8rgwyViYyzGdS502zXx4Tny1VdUdVXddOPwh8BdhtXLUZ+2zpUmDbDVjdU76NjQ/8R3Wq6lHgfmCHAdedq6Z7bK+nSfNjnphkNMkXkhy+KRo4Swbtl1e0p6EvSLLHNNedawY+rvay+V7A5T2z5+tYGcREfTdfx8p0jf9cKeCSJCuSLJulNs2m5yX57yQXJdmvnedYAZJsRRM8/rln9rwfL2lu0XoWcM24RTP22TJX/zSV+kjyKmAE+NWe2Qur6vYkTwMuT/Klqrpldlq42V0InFtVDyf5HZqzsy+a5TZ1xVLggqr6Qc+8x/JY0QSSvJAmsD2/Z/bz27HyFODSJF9tz8A8FlxH815Zl+SlwL8A+8xym7rkZcBVVdV7Nm5ej5ckW9ME1DdV1QObaj9dOsN2O7BHT3n3dl7fOkkWANsA9w647lw10LElORA4ATisqh4em19Vt7f/rgKuoPkNYD6Ysl+q6t6evvgo8JxB152jpnNcSxl3yWIej5VBTNR383WsDCTJM2jeO0uq6t6x+T1j5S7gU8yfW1CmVFUPVNW6dvrfgccn2ZHH+FjpMdlny7wbL0keTxPWzqmqT/apMnOfLbN9017PjXkLaG6624v1N2zuN67OsWz40MEn2un92PChg1XMn4cOBumXZ9Hc7LrPuPnbAVu20zsCNzNPboIdsF926Zl+OfCFdnp74Ott/2zXTm8/28e0OfqkrfczNDcB57EwVnqOcRET30j+a2x4Y/B/zuexMmCf7ElzP/AvjZv/JODJPdNXA4fM9rFsxn556th7hyZ4fLMdNwO9/+byz2T90i7fhuY+tyc9FsZL+7r/A3DaJHVm7LOlM5dEq+rRJMcBF9M8bXNGVd2Y5ERgtKqWAx8DzkqykmZQLG3XvTHJJ4AvA48Cx9aGl3rmrAH75a+ArYF/ap7B4JtVdRjws8DfJfkhzdnUU6rqy7NyIDNswH75/SSH0YyJNTRPjVJVa5KcBFzbbu7E2vD0/Zw0YJ9A8745r9pPjda8HSsASc6lebpvxyS3AW8HHg9QVR8G/p3maa6VwEPA0e2yeTlWYKA+eRvNPcIfaj9XHq2qEWBn4FPtvAXAP1bVZzb7AWwiA/TLEcDvJnkU+C6wtH0v9X3/zcIhbBID9As0vxhfUlXf6Vl1Po+XXwZeDXwpyfXtvLfS/LIz458t/mkqSZKkjuvSPWySJEnqw8AmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeq4/wGMoExHZjwsOQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "4wzFGvUrb5S7" }, "source": [ "## Exact Solution\n", "The initial value problem has the exact solution\n", "\\begin{equation}y(t)=2e^{-t}+t-1.\\end{equation}\n", "The figure below plots the exact solution." ] }, { "cell_type": "code", "metadata": { "id": "yNBE4FYob5S8", "colab": { "base_uri": "https://localhost:8080/", "height": 312 }, "outputId": "e79dcb6b-cfae-4829-c259-d006b3f3b243" }, "source": [ "IC=1 # Intial condtion\n", "y=(IC+1)*np.exp(-t)+t-1\n", "fig = plt.figure(figsize=(6,4))\n", "plt.plot(t,y,'o-',color='black')\n", "plt.title('Exact Solution ')\n", "plt.xlabel('time')" ], "execution_count": 5, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "Text(0.5, 0, 'time')" ] }, "metadata": {}, "execution_count": 5 }, { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEWCAYAAAB2X2wCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXhU9dnG8e8DsguCLBajJGJ5iyiKktpKRcAVkIIbAgaEiCxWVKTghoqVRpGyK4tBwYoREnGFIpRdXyhIBMFdURbDIggiFCQQ8nv/mInvGAJJYGbOLPfnuuZi5pwzc+6cDM/88pwz55hzDhERiX5lvA4gIiLBoYIuIhIjVNBFRGKECrqISIxQQRcRiREq6CIiMUIFXeQEmNlGM7v6BJ9bz8z+a2Zlg51L4psKuoSFvwD+7C9kBbfnQri+lmaWU8wyZ5nZ62b2g5n9ZGafmFmPEGT5VfF3zm12zp3qnDsS7HVJfDvF6wASV/7snFvgdYgA04C1QCKQCzQGfuNpIpGToBG6eM7MJprZ6wGPnzGzheZTw8xmm9lOM/vRf/+sgGVPN7OpZrbVP/8tM6sCvAucGfDXwJlFrPr3wEvOuf3OuTzn3Brn3LsBr93ezD41sz1mtsTMzjtG/pfM7O8Bj3/568DMpgH1gFn+HA+YWZKZOTM7xb/MmWb2jpntNrP1ZtYr4LWeMLMsM3vZzPb58ySf6LaW2KaCLpHgr0BjM+thZs2BnkB35zsvRRlgKr5RdD3gZyCwVTMNqAycD9QBRjvn9gNtgK3+1sapzrmtRax3BTDezDqbWb3AGWb2P8B0oD9QG5iDryiXL80P5pzrBmzG99fJqc654UUsNgPIAc4EbgGeMrMrA+a39y9THXin0M8v8gsVdAmnt/yj3YJbLwDn3AGgGzAKeAW4xzmX45+3yzn3unPugHNuH5AGtAAws7r4Cndf59yPzrnDzrmlpcjTEXgfeAzYYGYfmdnv/fM6Af9yzs13zh0GRgCVgGYnuQ1+xczOBv4EPOicO+ic+wh4Abg9YLH/dc7N8ffcpwEXBTODxA4VdAmnG5xz1QNukwtmOOdWAt8CBmQVTDezymb2vJltMrO9wHtAdf8RImcDu51zP55IGP+HwEPOufOBM4CP8H3oGL7R8qaAZfOB74CEE1nXcZyJ72fYFzBtU6H1bA+4fwCoWNCuEQmkgi4RwczuBioAW4EHAmb9Ffgd8AfnXDXgioKn4Cuwp5tZ9SJeslSnEXXO/YBvFH4mcLo/R2JAPsP3AbKliKfvx9f2KVB4x+rxsmzF9zNUDZhW7xjrETkuFXTxnL9f/XegK77WywNm1sQ/uyq+vvkeMzsdGFLwPOfcNnw7Pyf4d56WM7OCgv89UNPMTjvOep8xswvM7BR/Qb0LWO+c24Xvr4TrzewqMyuH74MlF1hexEt9BLT176D9Db6+e6DvgfpFZXDOfed/zafNrKKZXYhvH8Irx8otciwq6BJOBUd6FNze9LcOXgGecc6tdc59DTwCTDOzCsAYfL3rH/DtxJxb6DW7AYeBL4Ad+Iupc+4LfDs1v/X364s6yqUy8CawB1+7JxHfDkicc1/i+4B51r/uP+PbsXmoiNcpOPxxI/BvILPQ/KeBR/05Bhbx/C5AEr7R+pvAkAg7vFOihOkCFyIisUEjdBGRGKGCLiISI1TQRURihAq6iEiM8OzLCbVq1XJJSUlerV5EJCp9+OGHPzjnahc1z7OCnpSURHZ2tlerFxGJSma26Vjz1HIREYkRKugiIjFCBV1EJEaooIuIxAgVdBGRGKGCLiISJhkZGSQlJVGmTBmSkpLIyMgI6uvrJPkiImGQkZFB7969OXDgAACbNm2id+/eAKSkpARlHRqhi4iEweDBg38p5gUOHDjA4MGDg7YOFXQRkTDYvHlzqaafCBV0EZEwSEgo+nK09erVC9o6VNBFRMLg3HPPPWpa5cqVSUtLC9o6VNBFREJs7ty5LF26lPbt25OYmIiZkZiYSHp6etB2iIKOchERCam9e/fSq1cvzjvvPDIzM6lYsWLI1qWCLiISQoMGDWLr1q0sX748pMUc1HIREQmZBQsWkJ6ezl//+lf+8Ic/hHx95pwL+UqKkpyc7HQ+dBGJVfv27aNx48ZUqFCBjz76iEqVKgXldc3sQ+dcclHz1HIREQmBBx98kM2bN/O///u/QSvmxVHLRUQkyBYvXszEiRPp378/zZo1C9t61XIREQmi/fv307hxY8qWLcvatWupXLlyUF9fLRcRkTB5+OGH2bhxI0uXLg16MS+OWi4iIkHy/vvv8+yzz9KvXz+aN28e9vWr5SIiEgQHDhzgoosu4siRI3z88cdUqVIlJOtRy0VEJMQeffRR1q9fz6JFi0JWzIujlouIyElavnw5Y8aM4a677qJVq1ae5VBBFxE5CT///DOpqanUq1ePZ555xtMsarmIiJyEIUOG8NVXXzF//nyqVq3qaRaN0EVETtDKlSsZOXIkvXr14uqrr/Y6jgq6iMiJOHjwIKmpqSQkJDBixAiv4wBquYiInJAnn3ySzz//nLlz51KtWjWv4wAaoYuIlFp2djbDhw/njjvu4LrrrvM6zi+KLehmNsXMdpjZJ8eYn2Jm68zsYzNbbmYXBT+miEhkyM3NJTU1lTPOOIORI0d6HedXSjJCfwlofZz5G4AWzrnGwFAgPQi5REQiUlpaGp988gnp6elUr17d6zi/UmwP3Tn3npklHWf+8oCHK4CzTj6WiEjkWbNmDU899RS33347119/vddxjhLsHnpP4N1jzTSz3maWbWbZO3fuDPKqRURC59ChQ6SmplK7dm1Gjx7tdZwiBe0oFzNrha+gX36sZZxz6fhbMsnJyd6cFUxE5AQMGzaMtWvX8tZbb3H66ad7HadIQSnoZnYh8ALQxjm3KxivKSISKdatW8fQoUO57bbb6NChg9dxjumkWy5mVg94A+jmnPvq5COJiESOw4cPk5qayumnn864ceO8jnNcxY7QzWw60BKoZWY5wBCgHIBzbhLwOFATmGBmAHnHOleviEi0GT58OKtXr+b111+nZs2aXsc5Ll3gQkTkGD799FMuueQSbrjhBjIzM72OAxz/Ahf6pqiISBHy8vJITU2lWrVqPPfcc17HKRGdy0VEpAgjR45k1apVZGZmUrt2ba/jlIhG6CIihXzxxRcMGTKEm266iY4dO3odp8RU0EVEAhw5coTU1FSqVKnChAkT8B/sERXUchERCTBmzBhWrFhBRkYGZ5xxhtdxSkUjdBERv6+++opHH32U9u3b06VLF6/jlJoKuogIvlbLHXfcQcWKFZk0aVJUtVoKqOUiIgI899xzLFu2jH/+85/UrVvX6zgnRCN0EYl769ev5+GHH6Zt27Z069bN6zgnTAVdROJafn4+PXv2pHz58qSnp0dlq6WAWi4iEtcmTpzIe++9x4svvkhCQoLXcU6KRugiErc2bNjAgw8+yHXXXUdqaqrXcU6aCrqIxCXnHHfeeSdlypSJ+lZLAbVcRCQupaens2jRIp5//nnq1avndZyg0AhdROLOpk2bGDhwIFdffTW9evXyOk7QqKCLSFxxztGrVy+cc0yePDkmWi0F1HIRkbgyZcoU5s+fz/jx40lKSvI6TlBphC4icSMnJ4cBAwbQsmVL+vbt63WcoFNBF5G44Jyjd+/e5OXl8eKLL1KmTOyVP7VcRCQuvPzyy7z77ruMHTuW+vXrex0nJGLvI0pEpJCtW7fSv39/mjdvTr9+/byOEzIq6CIS05xz9OnTh4MHD8Zsq6WAWi4iEtMyMjKYPXs2o0aNokGDBl7HCanY/agSkbi3fft27r33Xi677DLuvfder+OEnAq6iMQk5xx33XUXBw4cYMqUKZQtW9brSCGnlouIxKTMzEzeeusthg8fTsOGDb2OExYaoYtIzNmxYwf9+vXj0ksvZcCAAV7HCZtiC7qZTTGzHWb2yTHmNzSz/5hZrpkNDH5EEZHSufvuu9m3bx9Tp06Ni1ZLgZKM0F8CWh9n/m7gXmBEMAKJiJyMmTNnMnPmTJ544gkaNWrkdZywKragO+few1e0jzV/h3NuFXA4mMFERErrhx9+4C9/+QtNmzZl0KBBXscJu7DuFDWz3kBvIGZOKC8ikeOee+5hz549LFy4kFNOib9jPsK6U9Q5l+6cS3bOJdeuXTucqxaRGPfWW28xY8YMHnvsMRo3bux1HE/oKBcRiXq7d++mb9++NGnShIceesjrOJ6Jv79JRCTm3HfffezatYu5c+dSrlw5r+N4ptiCbmbTgZZALTPLAYYA5QCcc5PM7DdANlANyDez/kAj59zekKUWEfGbNWsWr7zyCo8//jhNmjTxOo6nzDnnyYqTk5Nddna2J+sWkdjw448/cv7551OrVi2ys7MpX76815FCzsw+dM4lFzVPLRcRiVoDBgxgx44dzJ49Oy6KeXG0U1REotK7777LSy+9xIMPPsgll1zidZyIoJaLiESdn376ifPPP5/TTjuN1atXU6FCBa8jhY1aLiISUwYOHMi2bdt4880346qYF0ctFxGJKv/+97954YUXGDhwIL///e+9jhNR1HIRkaixb98+LrjgAipXrsyaNWuoWLGi15HCTi0XEYkJDzzwAN999x3Lli2Ly2JeHLVcRCQqLFq0iEmTJjFgwAAuu+wyr+NEJLVcRCTi/fe//6Vx48aUK1eOtWvXUqlSJa8jeUYtFxGJag8//DCbNm3ivffei+tiXhy1XEQkoi1dupTnnnuOe++9l8svv9zrOBFNLRcRiVj79+/noosuwjnHunXrqFKliteRPKeWi4hEpcGDB/PNN9+wePFiFfMSUMtFRCLSsmXLGDduHHfffTctW7b0Ok5UUEEXkYjz888/k5qaSmJiIsOGDfM6TtRQy0VEIs5jjz3G119/zcKFCzn11FO9jhM1NEIXkYiyYsUKRo8eTZ8+fbjyyiu9jhNVVNBFJGIcPHiQ1NRUEhISGD58uNdxoo5aLiISMZ544gm++OIL5s2bR7Vq1byOE3U0QheRiLBq1Sr+8Y9/0LNnT6699lqv40QlFXQR8Vxubi49evSgbt26jBw50us4UUstFxHx3NChQ/nss8/417/+xWmnneZ1nKilEbqIeGr16tUMGzaM7t2707ZtW6/jRDUVdBHxzKFDh+jRowd16tRh9OjRXseJemq5iIhnnnrqKT7++GPeeecdatSo4XWcqKcRuoh4Yu3ataSlpZGSksKf//xnr+PEhGILuplNMbMdZvbJMeabmY0zs/Vmts7MLgl+TJ+MjAySkpIoU6YMSUlJZGRkhGpVIhJChw8fpkePHtSsWZOxY8d6HSdmlGSE/hLQ+jjz2wAN/LfewMSTj3W0jIwMevfuzaZNm3DOsWnTJnr37q2iLhKFnnnmGT766CMmTpxIzZo1vY4TM0p0gQszSwJmO+cuKGLe88AS59x0/+MvgZbOuW3He83SXuAiKSmJTZs2HTU9MTGRjRs3lvh1RMRbH3/8MU2bNuXmm29m+vTpXseJOse7wEUweugJwHcBj3P804oK0tvMss0se+fOnaVayebNm0s1XUQiT15eHqmpqVSvXp1nn33W6zgxJ6w7RZ1z6c65ZOdccu3atUv13Hr16pVquohEnhEjRvDhhx8yfvx4atWq5XWcmBOMgr4FODvg8Vn+aUGVlpZG5cqVj5p+3333BXtVIhICn332GUOGDOGWW26hY8eOXseJScEo6O8At/uPdvkj8FNx/fMTkZKSQnp6OomJiZgZCQkJVKpUiVdeeYXc3Nxgr05EgujIkSPccccdVK1alfHjx3sdJ2aV5LDF6cB/gN+ZWY6Z9TSzvmbW17/IHOBbYD0wGfhLqMKmpKSwceNG8vPzycnJYfr06axevZqBAweGapUiEgSjR49m5cqVPPvss9SpU8frODGrREe5hEJpj3I5lgEDBjB69GhmzpzJzTffHIRkIhJMX375JRdddBFt2rThjTfewMy8jhTVjneUS9QX9EOHDtG8eXO++OIL1qxZQ/369YOQTkSC4ciRI1xxxRV8/vnnfPbZZ/zmN7/xOlLUC/Vhi54qX748M2bMoEyZMnTq1En9dJEIMm7cOJYvX864ceNUzMMg6gs6wDnnnMPUqVPJzs7mgQce8DqOiADr169n8ODBtGvXjpSUFK/jxIWYKOgAN9xwA/fddx/jxo3jzTff9DqOSFzLz8/njjvuoHz58kyaNEl98zCJmYIOMHz4cJKTk0lNTWXDhg1exxGJW+PHj+f9999nzJgxJCQU+cVxCYGYKujly5cnMzMTgE6dOnHo0CGPE4nEn2+//ZaHHnqINm3a0L17d6/jxJWYKugA9evXZ8qUKaxatYoHH3zQ6zgicSU/P5+ePXtStmxZnn/+ebVawizmCjrATTfdxD333MOYMWN4++23vY4jEjeef/55lixZwqhRozj77LOLf4IEVdQfh34subm5/OlPf+Kbb75hzZo1JCUlhWxdIgIbN26kcePGXHbZZcybN0+j8xCJ6ePQj6VChQpkZmaSn59P586d1U8XCSHnHL169QJg8uTJKuYeidmCDnDuuefy4osvsnLlSh5++GGv44jEnMDLQi5YsICOHTuSmJjoday4FdMFHeCWW27h7rvvZtSoUbzzzjtexxGJGYGXhSyQmZmpy0J6KGZ76IEOHjxIs2bN2LhxI2vWrNEIQiQIdFlIb8RlDz1QxYoVycrKIi8vj86dO3P48GGvI4lEPV0WMvLERUEH+O1vf8sLL7zAihUreOSRR7yOIxLV8vLyqFSpUpHzdFlI78RNQQe49dZbueuuuxgxYgSzZ8/2Oo5IVDp8+DBdunThwIEDlCtX7lfzKleuTFpamkfJJK4KOsCoUaNo0qQJ3bt357vvvvM6jkhUOXToEJ06dWLmzJmMGjWKqVOn/nJZyMTERNLT03VmRQ/FxU7Rwr7++msuueQSLrzwQpYsWXLUKENEjpabm0vHjh2ZNWsW48aN45577vE6UlyK+52ihTVo0IDJkyezfPlyHn30Ua/jiES8gwcPcuONNzJr1iwmTpyoYh6h4rKgA3Tu3Jk+ffowfPhw5syZ43UckYh14MAB2rdvz9y5c5k8eTJ9+/Yt/kniibgt6OC7EvmFF17I7bffTk5OjtdxRCLO/v37adeuHQsWLGDq1KnceeedXkeS44jrgl6pUiWysrLIzc2lc+fO5OXleR1JJGLs27ePtm3bsnTpUqZNm6Zzm0eBuC7oAL/73e94/vnnWbZsGY899pjXcUQiwt69e2ndujXLli3j1Vdf1ZErUSLuCzrAbbfdRq9evRg2bBhz5871Oo6Ip/bs2cO1117LBx98QGZmJp06dfI6kpSQCrrf2LFjady4Md26dWPLli1exxHxxO7du7nmmmtYvXo1M2fO5Oabb/Y6kpSCCrpfQT/9559/pkuXLuqnS9zZtWsXV111FevWreONN96gQ4cOXkeSUipRQTez1mb2pZmtN7OHipifaGYLzWydmS0xs7OCHzX0GjZsyKRJk3j//fcZMmSI13FEwmbnzp20atWKzz//nLfffpt27dp5HUlOQLEF3czKAuOBNkAjoIuZNSq02AjgZefchcCTwNPBDhouXbt2pWfPnjz99NPMmzfP6zgiIff999/TqlUr1q9fz+zZs2ndurXXkeQElWSEfimw3jn3rXPuEDADKPy3WCNgkf/+4iLmR5Vx48Zx/vnn061bN7Zu3ep1HJGQ2bZtGy1btmTDhg3MmTOHq6++2utIchJKUtATgMCzWOX4pwVaC9zkv38jUNXMap58PG9UrlyZrKws9u/fz2233aZ+usSknJwcWrRoQU5ODnPnzqVly5ZeR5KTFKydogOBFma2BmgBbAGOFF7IzHqbWbaZZe/cuTNIqw6N8847j4kTJ7J06VL+9re/eR1HJKg2b95MixYt2L59O/PmzaN58+ZeR5IgKElB3wKcHfD4LP+0XzjntjrnbnLOXQwM9k/bU/iFnHPpzrlk51xy7dq1TyJ2eNx+++2kpqaSlpbG/PnzvY4jEhQbN26kRYsW7Nq1iwULFtCsWTOvI0mQlKSgrwIamNk5ZlYe6Az86mrLZlbLzApe62FgSnBjeufZZ5/lvPPOIyUlhW3btnkdR+SkfPPNN1xxxRX89NNPLFy4kEsvvdTrSBJExRZ051we0A+YB3wOZDnnPjWzJ82svX+xlsCXZvYVcAYQM5csqVKlCq+99tov/fQjR47qJIlEha+++ooWLVpw4MABFi1aRNOmTb2OJEEWlxe4OBEvvfQSqampPP744+qpS9T54osvuPLKK8nLy2PhwoU0btzY60hygnSBiyDo0aMH3bt3Z+jQoSxcuNDrOCIl9umnn9KiRQvy8/NZsmSJinkMU0EvhfHjx9OwYUNSUlLYvn2713FEirVu3TpatmxJ2bJlWbJkCY0aFf5OoMQSFfRSqFKlCllZWezdu5eUlBT10yWirVmzhlatWlGxYkWWLl1Kw4YNvY4kIaaCXkoXXHABzz33HIsWLeLvf/+713FEipSdnc2VV17JqaeeytKlS2nQoIHXkSQMVNBPQGpqKt26deNvf/sbixYtKv4JImG0cuVKrr76amrUqMHSpUupX7++15EkTFTQT4CZMWHCBH73u9+RkpLC999/73UkEQCWLVvGNddcQ61atVi6dClJSUleR5IwUkE/QaeeeipZWVns2bOHrl27qp8unnvvvfe47rrrqFu3LkuXLuXss88u/kkSU1TQT0Ljxo159tlnWbBgAU899ZTXcSSOLVq0iDZt2lCvXj2WLFlCQkLh8+dJPFBBP0k9e/YkJSWFJ554giVLlngdR+LQv//9b66//nrq16/P4sWLqVu3rteRxCMq6CfJzJg4cSK//e1vue2229ixY4fXkSSOzJkzh/bt2/M///M/LFq0iDPOOMPrSOIhFfQgqFq1Kq+99ho//vgjXbt2JT8/3+tIEgdmzZrFjTfeSKNGjVi0aBHRcAZTCS0V9CC58MILGTt2LPPnz+fpp6P2CnwSJd58801uuukmLrroIhYuXEjNmlF7PRkJIhX0IOrVqxddunTh8ccf57333vM6jsSo1157jY4dO5KcnMz8+fOpUaOG15EkQqigB5GZ8fzzz3PuuefSpUsXIv2qTBJ9pk+fTpcuXfjjH//IvHnzOO2007yOJBFEBT3IqlatSlZWFrt27aJbt27qp0vQTJs2ja5du3L55Zczd+5cqlWr5nUkiTAq6CHQpEkTxowZw7x583jmmWe8jiMxYMqUKXTv3p1WrVoxZ84cTj31VK8jSQRSQQ+RPn360KlTJx577DHef/99r+NIFEtPT6dnz55ce+21zJo1i8qVK3sdSSKUCnqImBnp6emcc845dOnShR9++MHrSBKFxo8fT58+fbj++ut56623qFSpkteRJIKpoIdQtWrVyMrKYufOndx+++3qp0upjB07ln79+tGhQwdef/11Klas6HUkiXAq6CF28cUXM3r0aN59913+8Y9/eB1HosSIESPo378/N910E1lZWVSoUMHrSBIFVNDD4K677qJjx44MHjyYZcuWeR1HItzTTz/NoEGDuPXWW5kxYwbly5f3OpJECRX0MDAzJk+eTGJiIp07d2bXrl1eR5II9eSTT/LII49w2223kZGRQbly5byOJFFEBT1MTjvtNLKystixYwfdu3dXP11+xTnH448/zpAhQ+jevTsvv/wyp5xyitexJMqooIdR06ZNGTlyJP/6178YOXKk13EkQjjneOSRRxg6dCg9e/ZkypQplC1b1utYEoVU0MPs7rvv5uabb+bhhx9m+fLlXscRjznnGDRoEMOGDaNv376kp6dTpoz+W8qJ0TsnzMyMF198kXr16qmfHuecc9x///2MHDmSfv36MWHCBBVzOSl693igoJ++fft2evTogXPO60gSZvn5+fTr14+xY8dy//33M27cOMzM61gS5UpU0M2stZl9aWbrzeyhIubXM7PFZrbGzNaZWdvgR40tycnJjBgxgtmzZzNq1Civ40gY5efn07dvXyZMmMADDzzAyJEjVcwlOJxzx70BZYFvgPpAeWAt0KjQMunAXf77jYCNxb1u06ZNXbzLz893N954ozvllFPcf/7zH6/jSBjk5eW51NRUB7jBgwe7/Px8ryNJlAGy3THqaklG6JcC651z3zrnDgEzgA6FPxeAgnN5ngZsPYnPmLhhZkyZMoWzzjqLTp06sXv3bq8jSQgdOXKE1NRUpk6dyhNPPMHQoUM1MpegKklBTwC+C3ic458W6Amgq5nlAHOAe4p6ITPrbWbZZpatiz/4VK9enczMTLZt20Zqaqr66TEqLy+Prl27Mm3aNP7+978zZMgQFXMJumDtFO0CvOScOwtoC0wzs6Ne2zmX7pxLds4l64K2/+/SSy9l+PDhvPPOO4wZM8brOBJkhw8fpkuXLsyYMYNnnnmGwYMHex1JYlRJCvoW4OyAx2f5pwXqCWQBOOf+A1QEagUjYLy477776NChAw8++CAffPCB13EkSA4dOkSnTp2YOXMmo0aN4oEHHvA6ksSwkhT0VUADMzvHzMoDnYF3Ci2zGbgKwMzOw1fQ1VMpBTNj6tSpnHnmmdx66638+OOPXkeSk5Sbm8stt9zCm2++ybhx47j//vu9jiQxrtiC7pzLA/oB84DPgSzn3Kdm9qSZtfcv9legl5mtBaYDPZyawaVWo0YNMjMz2bJlC3fccYf66VHs4MGD3HjjjcyaNYuJEydyzz1F7lYSCa5jHf4S6psOWzy2UaNGOcCNGTPG6yhyAvbv3++uueYaZ2Zu8uTJXseRGMNJHrYoYda/f3/at2/PoEGDWLVqlddxpBT2799Pu3btWLBgAVOnTuXOO+/0OpLEERX0CFTQT69bty633nore/bs8TqSlMC+ffto27YtS5cuZdq0aXTv3t3rSBJnVNAj1Omnn05mZiY5OTn07NlT/fQIt3fvXlq3bs2yZct49dVXSUlJ8TqSxCEV9Aj2xz/+kWHDhvHGG2/w3HPPeR1HjmHPnj1ce+21fPDBB2RmZtKpUyevI0mcUkGPcAMGDPoA6UEAAAyGSURBVKBdu3YMHDiQ7Oxsr+NIIbt37+aaa65h9erVzJw5k5tvvtnrSBLHVNAjnJnx0ksvccYZZ9CpUyd++uknryOJ365du7jqqqtYt24db7zxBh06FD7FkUh4qaBHgZo1azJjxgw2bdrEnXfeqX56BNi5cyetWrXi888/5+2336Zdu3ZeRxJRQY8WzZo14+mnn2bmzJlMmDDB6zhxbfv27bRs2ZL169cze/ZsWrdu7XUkEUAFPar89a9/pW3btgwYMIDVq1d7HScubd26lZYtW7Jx40bmzJnD1Vdf7XUkkV+ooEeRMmXK8M9//pM6depw6623snfvXq8jxZWcnBxatmzJli1bmDt3Li1btvQ6ksivqKBHmVq1ajFjxgw2btxIr1691E8Pk02bNtGiRQu2b9/OvHnzaN68udeRRI6igh6F/vSnP5GWlkZWVhaTJk3yOk7M27BhAy1atGDXrl0sWLCAZs2aeR1JpEgq6FFq0KBBtGnThv79+7NmzRqv48Ss9evX06JFC/bu3cvChQu59NJLvY4kckwq6FGqoJ9eu3Zt9dND5KuvvqJFixYcOHCARYsW0bRpU68jiRyXCnoUq127NtOnT2fDhg307t1b/fQgyMjIICkpiTJlynDeeeexb98+Fi9eTJMmTbyOJlIsFfQo17x5c4YOHUpmZibp6elex4lqGRkZ9O7dm02bNuGcIz8/n8OHD7Nu3Tqvo4mUiHk1qktOTnY6N0lw5Ofn07ZtW5YsWcKKFSs0mizGwYMH2bJly1G39PR0Dhw4cNTyiYmJbNy4MfxBRYpgZh8655KLnKeCHht27NjBxRdfTJUqVfjwww+pWrWq15HCzjnHrl27finQOTk5RRbu3bt3H/XcKlWqsH///iJf18zIz88PdXyREjleQT8l3GEkNOrUqcP06dNp1aoVffr0ISMjAzPzOlbQ5ObmsnXr1iILdMFt69at5Obm/up5ZkadOnVISEggKSmJyy+/nISEhKNu1apV45xzzmHTpk1HrbtevXrh+jFFTooKegy54oorePLJJ3n00Udp1aoVvXr18jpSsZxz/Pjjj8cs0gWj7B9++OGo51aqVOmXgnzZZZcVWajr1q1LuXLlSpQlLS2N3r17/6rtUrlyZdLS0oL284qEklouMSY/P5/WrVvz/vvvs3LlSi688ELPshw+fJht27Yds/VRcDt48OBRz61du3aRBTohIYGzzjqLhIQEqlevHvS/QjIyMhg8eDCbN2+mXr16pKWl6epDElHUQ48zO3bsoEmTJlStWpXs7Oyg99Odc/z000/HLdJbtmxhx44dRx1KWaFCBc4888yjinPhUXWFChWCmlkkVqiHHmfq1KnDq6++ylVXXUXr1q3Jycnhu+++K9GIMy8vj+3btx+z9VFwK+pokJo1a/5SlC+55JIiR9c1a9aMqd6+SCTRCD2G3XLLLbz++uu/mlaxYkX69+9Pw4YNixxVf//990cd0VGuXLlfjaqLup155plUqlQpnD+eSFzSCD1OrVq16qhpBw8eZNiwYb88rl69+i9FuXHjxkW2QWrVqkWZMvoOmkikU0GPYd99912R082ML7/8koSEBCpXrhzmVCISKiUadplZazP70szWm9lDRcwfbWYf+W9fmdme4EeV0jrW8dP16tWjQYMGKuYiMabYgm5mZYHxQBugEdDFzBoFLuOcu98518Q51wR4FngjFGGldNLS0o4q2jquWiR2lWSEfimw3jn3rXPuEDAD6HCc5bsA04MRTk5OSkoK6enpJCYmYmYkJiaSnp6u46pFYlRJeugJQGAzNgf4Q1ELmlkicA6w6OSjSTCkpKSogIvEiWAfutAZmOmcO1LUTDPrbWbZZpa9c+fOIK9aRCS+laSgbwHODnh8ln9aUTpznHaLcy7dOZfsnEuuXbt2yVOKiEixSlLQVwENzOwcMyuPr2i/U3ghM2sI1AD+E9yIIiJSEsUWdOdcHtAPmAd8DmQ55z41syfNrH3Aop2BGU7XQRMR8USJvljknJsDzCk07fFCj58IXiwRESktz87lYmY7gaOvJlAytYCjT5DtvUjNBZGbTblKR7lKJxZzJTrnitwJ6VlBPxlmln2sk9N4KVJzQeRmU67SUa7SibdcOuOSiEiMUEEXEYkR0VrQ070OcAyRmgsiN5tylY5ylU5c5YrKHrqIiBwtWkfoIiJSiAq6iEiMiLiCXoKLaVQws0z//JVmlhQw72H/9C/N7Low5xpgZp+Z2TozW+g/82TBvCMBFwA56rQJIc7Vw8x2Bqz/zoB53c3sa/+te5hzHfOiKCHeXlPMbIeZfXKM+WZm4/y515nZJQHzQrm9isuV4s/zsZktN7OLAuZt9E//yMyCeqHeEuRqaWY/Bfy+Hg+Yd9z3QIhzDQrI9In/PXW6f15ItpeZnW1mi/114FMzu6+IZUL7/nLORcwNKAt8A9QHygNrgUaFlvkLMMl/vzOQ6b/fyL98BXyn8P0GKBvGXK2Ayv77dxXk8j/+r4fbqwfwXBHPPR341v9vDf/9GuHKVWj5e4Apod5e/te+ArgE+OQY89sC7wIG/BFYGertVcJczQrWh+9iMysD5m0Eanm0vVoCs0/2PRDsXIWW/TOwKNTbC6gLXOK/XxX4qoj/jyF9f0XaCL0kF9PoAPzTf38mcJWZmX/6DOdcrnNuA7De/3phyeWcW+ycO+B/uALfWSlDrbQXHwl0HTDfObfbOfcjMB9o7VGusF0UxTn3HrD7OIt0AF52PiuA6mZWl9Bur2JzOeeW+9cL4Xt/lWR7HcvJvDeDnSss7y/n3Dbn3Gr//X34zn2VUGixkL6/Iq2gF3UxjcIb5JdlnO/EYT8BNUv43FDmCtQT36dwgYrmOw/8CjO7IUiZSpPrZv+fdzPNrOBUyBGxvazoi6KEanuVxLGyh3J7lVbh95cD/m1mH5pZbw/yXGZma83sXTM73z8tIraXmVXGVxhfD5gc8u1lvlbwxcDKQrNC+v4q0cm5pOTMrCuQDLQImJzonNtiZvWBRWb2sXPumzBFmgVMd87lmlkffH/dXBmmdZdEURdF8XJ7RTQza4WvoF8eMPly//aqA8w3sy/8I9hwWI3v9/VfM2sLvAU0CNO6S+LPwDLnXOBoPqTby8xOxfcB0t85tzdYr1sSkTZCL8nFNH5ZxsxOAU4DdpXwuaHMhZldDQwG2jvncgumO+e2+P/9FliC75M7LLmcc7sCsrwANC3pc0OZK8BRF0UJ4fYqiWNlD+X2KhEzuxDf77CDc25XwfSA7bUDeJPgtRqL5Zzb65z7r//+HKCcmdUiAraX3/HeX0HfXmZWDl8xz3DOvVHEIqF9fwV7x8BJ7lQ4Bd/OgHP4/x0p5xda5m5+vVM0y3//fH69U/RbgrdTtCS5Lsa3E6hBoek1gAr++7WArwnSzqES5qobcP9GYIX7/50wG/z5avjvnx6uXP7lGuLbQWXh2F4B60ji2Dv5rufXO60+CPX2KmGuevj2CzUrNL0KUDXg/nKgdRhz/abg94evMG72b7sSvQdClcs//zR8ffYq4dhe/p/7ZWDMcZYJ6fsraBs3iL+ktvj2Dn8DDPZPexLfqBegIvCa/839AVA/4LmD/c/7EmgT5lwLgO+Bj/y3d/zTmwEf+9/QHwM9w5zraeBT//oXAw0DnnuHfzuuB1LDmcv/+AlgWKHnhXp7TQe2AYfx9Sl7An2Bvv75Boz35/4YSA7T9iou1wvAjwHvr2z/9Pr+bbXW/3seHOZc/QLeXysI+MAp6j0Qrlz+ZXrgO1Ai8Hkh21742mAOWBfwe2obzveXvvovIhIjIq2HLiIiJ0gFXUQkRqigi4jECBV0EZEYoYIuIhIjVNAlbphZdTP7i//+mWY20+tMIsGkwxYlbvjPrzHbOXeBx1FEQkLncpF4Mgw418w+wvcN1POccxeYWQ/gBnzfHGwAjMD37cZuQC7Q1jm328zOxfelkNrAAaCXc+6L8P8YIkVTy0XiyUPAN865JsCgQvMuAG4Cfg+kAQeccxcD/wFu9y+TDtzjnGsKDAQmhCW1SAlphC7is9j5zmG9z8x+wneWSvB9PftC/xn0mgGv+U6/D/jOGyQSMVTQRXxyA+7nBzzOx/f/pAywxz+6F4lIarlIPNmH79JgpeZ857XeYGYd4ZdrQ15UzNNEwkoFXeKG851DfJn/wsL/OIGXSAF6mlnBmfqCdkk1kWDQYYsiIjFCI3QRkRihgi4iEiNU0EVEYoQKuohIjFBBFxGJESroIiIxQgVdRCRG/B+yU8P/8f03ZgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "A6YrDJJnb5S-" }, "source": [ "## 2-step Adams Moulton\n", "\n", "The general 2-step Adams Moulton difference equation is\n", "\\begin{equation}w_{i+1} = w_{i} + \\frac{h}{12}(5f(t_{i+1},w_{i+1})+8f(t_{i},w_{i})-f(t_{i-1},w_{i-1})). \\end{equation}\n", "For the specific intial value problem the 2-step Adams Moiulton difference equation is\n", "\\begin{equation}w_{i+1} = w_{i} + \\frac{h}{12}(5(t_{i+1}-w_{i+1})+8(t_{i}-w_{i})-(t_{i-1}-w_{i-1})). \\end{equation}\n", "\n", "for $i=0$ the difference equation is:\n", "\\begin{equation}w_{1} = w_{0} + \\frac{h}{12}(5(t_{1}-w_{1})+8(t_{0}-w_{0})-(t_{-1}-w_{-1})).\\end{equation}\n", "\n", "this is not solvable as $w_{1}, \\ w_{-1}$ are unknown.\n", "for $i=1$ the difference equation is:\n", "\\begin{equation}w_{2} = w_{1} + \\frac{h}{12}(5(t_{2}-w_{2})+8(t_{1}-w_{1})-(t_{0}-w_{0})). \\end{equation}\n", "this is not solvable as $w_{1}$ and $w_{2}$ are unknown. $w_1$ can be approximated using a one step method. Here, as the exact solution is known,\n", "\\begin{equation}w_1=2e^{-t_1}+t_1-1.\\end{equation}\n", "As the intial value problem is linear the difference equation can be rearranged such that $w_2$ is on the right hand side:\n", "\\begin{equation}w_{2}+\\frac{5h}{12}w_{2} = w_{1} + \\frac{h}{12}(5(t_{2})+8f(t_{1}-w_{1})-(t_{0}-w_{0})), \\end{equation}\n", "\\begin{equation}w_{2} = \\frac{w_{1} + \\frac{h}{12}(5(t_{2})+8f(t_{1}-w_{1})-(t_{0}-w_{0}))}{1+\\frac{5h}{12}}. \\end{equation}\n" ] }, { "cell_type": "code", "metadata": { "id": "6FzTfyidb5S-" }, "source": [ "### Initial conditions\n", "w=np.zeros(len(t))\n", "w[0]=IC\n", "w[1]=y[1] # NEEDED FOR THE METHOD" ], "execution_count": 6, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "0cwXy_Enb5TA" }, "source": [ "### Loop" ] }, { "cell_type": "code", "metadata": { "id": "lNn9_F6Sb5TA" }, "source": [ "for k in range (1,N):\n", " w[k+1]=(w[k]+h/12.0*(5*t[k+1]+8*myfun_ty(t[k],w[k])-myfun_ty(t[k-1],w[k-1])))/(1+5*h/12) " ], "execution_count": 7, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "sZzEevGlb5TB" }, "source": [ "### Plotting solution" ] }, { "cell_type": "code", "metadata": { "id": "zTE-KOn_b5TB" }, "source": [ "def plotting(t,w,y):\n", " fig = plt.figure(figsize=(10,4))\n", " plt.plot(t,y, 'o-',color='black',label='Exact')\n", " plt.plot(t,w,'s:',color='blue',label='Adams-Moulton')\n", " plt.xlabel('time')\n", " plt.legend()\n", " plt.show " ], "execution_count": 8, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "GzHPj5j5b5TD" }, "source": [ "The plot below shows the exact solution (black) and the 2 step Adams-Moulton approximation (red) of the intial value problem" ] }, { "cell_type": "code", "metadata": { "id": "7thBBKleb5TE", "colab": { "base_uri": "https://localhost:8080/", "height": 283 }, "outputId": "36b380f6-3657-4ad8-e9d1-c8792890c132" }, "source": [ "plotting(t,w,y)" ], "execution_count": 9, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAEKCAYAAAA2KNBKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3yN5//H8deVCEJTiVHamm3tFQSt2qp0otU2qNIv1dboHjq1VIdqqaL4omgJQY3avvYMsfeqXaskRgmSXL8/cvjFTMhJ7pyc9/PxOA/nHuec95U7h4/rvu7rNtZaREREROT2+DgdQERERMSTqZgSERERSQEVUyIiIiIpoGJKREREJAVUTImIiIikgIopERERkRRIspgyxgw1xhw1xmy8wfZGxpj1xpi1xphIY0x198cUERERSZ9MUvNMGWNqAmeAEdbaMtfZfgfwr7XWGmPKAeHW2hKpklZEREQknUmyZ8pauxA4cZPtZ+z/V2TZAc0CKiIiIl4jkzvexBjTBPgGuAt4IjmvyZ07ty1cuLA7Pl5EREQkVa1ateofa22e621zSzFlrZ0ATHCdEuwGPHK9/Ywx7YB2AAULFiQyMtIdHy8iIiKSqowxe2+0za1X87lOCd5njMl9g+2DrLUh1tqQPHmuW9yJiIiIeJQUF1PGmAeMMcb1vCKQBTie0vcVERER8QRJnuYzxoQBtYHcxpgDQBfAD8BaOwB4FnjJGHMROAe8YJO6RFBEREQkg0iymLLWNkti+3fAd+4Ic/HiRQ4cOEBMTIw73k4ckDVrVvLnz4+fn5/TUURERNKEWwagu8uBAwcICAigcOHCuM4cigex1nL8+HEOHDhAkSJFnI4jIiKSJtLV7WRiYmLIlSuXCikPZYwhV65c6lkUEZFUly8fGHPtI1++tM+SroopQIWUh9PxExGRtHDkyK2tT03prpgSERER8SQqpq7i6+tLcHDw5ce3337rtvdeu3Yt06ZNc9v7iYiIiPM8upgaOXIkhQsXxsfHh8KFCzNy5MgUv6e/vz9r1669/OjcubMbkiZQMSUiIpLxeGwxNXLkSNq1a8fevXux1rJ3717atWvnloLqaidPnqR48eJs27YNgGbNmvHf//4XgNdff52QkBBKly5Nly5dLr9m5cqVVKtWjfLly1OlShVOnjzJ559/zpgxYwgODmbMmDFuzykiIuINTp8+7XSEK6SrqRESe+utt1i7du0Nty9fvpzz589fse7s2bO0adPmcqFzteDgYHr37n3Tzz137hzBwcGXlz/66CNeeOEF+vbtS+vWrXnzzTeJiorilVdeAaB79+7kzJmTuLg46tWrx/r16ylRogQvvPACY8aMoXLlypw6dYps2bLRtWtXIiMj6du3b3J/DCIiIuJy+DDMmBFN376PAFOAay/dy5HjHOCfprnSbTGVlKsLqaTWJ9el03xXq1+/PmPHjqVDhw6sW7fu8vrw8HAGDRpEbGwshw4dYvPmzRhjuPvuu6lcuTIAd955Z4oyiYiICLzxxinGj/chc+ZDTJmyiujoaD755BP27dtHwYIF6d69Oy1atEjzXOm2mEqqB6lw4cLs3XvtDZwLFSrE/Pnz3Z4nPj6eLVu2kC1bNqKiosifPz+7d++mZ8+erFy5kqCgIFq3bq05lkRERFLBmjVrmDfvWe64oyAzZ47nwQcfBHCkeLqax46Z6t69O9myZbtiXbZs2ejevXuqfF6vXr0oWbIko0aN4uWXX+bixYucOnWK7NmzkyNHDo4cOcL06dMBKF68OIcOHWLlypVAwrnd2NhYAgIC0t15XhERkfRs1iyoXfsINWvWI1u2OCIiBlwupNILjy2mWrRowaBBgyhUqBDGGAoVKsSgQYNSXKFeGjN16dG5c2e2bdvG4MGD+eGHH6hRowY1a9bkq6++onz58lSoUIESJUrQvHlzHn74YQAyZ87MmDFj6NSpE+XLl6d+/frExMRQp04dNm/erAHoIiIiyTR27AoWLDhM/vylWbp0KSVKlHA60jWMtdaRDw4JCbGRkZFXrNuyZQslS5Z0JI+4j46jiIik1KFDMHZsH958802qV6/Ln3+OJzAw0LE8xphV1tqQ623z2J4pERERyZi++85y331nefPNnjRp0oTZs6c6WkglJd0OQBcRERHvc/HiRSIiPiYmJpB27Z6gf/+++Pr6Oh3rptQzJSIiIo47cwYGDjxP48aNmTChJ19+GceAAf3TfSEF6pkSERGRdKBHj3/p1i0LxuxlwIABvPrqq05HSjYVUyIiIuKovXv3MmbMY/j55SI8/CsaN27sdKRbomJKREREHLFmDXTocJrduxsSE3OYOXMGUqNGDadj3TKNmbqOiRMnYoxh69at191eu3Ztrp7WIbXMnz8fYwyDBw++vG7t2rUYY+jZs+dtv+eTTz55+fnSpUvdklVERORWzJ+/noiIo1ibm0WLFnlkIQUeXEzlywfGXPvId+09D29ZWFgY1atXJywsLOVv5gZlypQhPDz88nJYWBjly5d3y3urmBIRkbR28CD88ccffPRRFYoWbcyKFSMpU6aM07Fum8cWU0eO3Nr65Dpz5gyLFy9myJAhjB49GkiYFT00NJSSJUvSpEkTzp07d3n/119/nZCQEEqXLk2XLl0ury9cuDAfffQRwcHBhISEsHr1aho0aMD999/PgAEDADh06BA1a9YkODiYMmXKsGjRoutmKlSoEDExMRw5cgRrLTNmzOCxxx67vH3t2rU8+OCDlCtXjiZNmhAVFQVc2YP2zz//ULhw4Sved8+ePQwYMIBevXoRHBzMokWL2LNnD3Xr1qVcuXLUq1ePffv2AdC6dWveeOMNqlWrxn333ce4ceNS9oMWERGvNG4cFCkSS9OmPalQoQJLlsynYMGCTsdKkXRdTNWuDcOGJTy/eDFh+fffk/faf/5J2P/PPxOWDx9O3usmTZpEw4YNKVasGLly5WLVqlX88ssvZMuWjS1btvDll1+yatWqy/t3796dyMhI1q9fz4IFC1i/fv3lbQULFmTt2rXUqFGD1q1bM27cOJYvX3656Bo1ahQNGjRg7dq1rFu3juDg4Bvmatq0KWPHjmXp0qVUrFiRLFmyXN720ksv8d1337F+/XrKli3Ll19+may2Fi5cmNdee4233377cs5OnTrRqlUr1q9fT4sWLXjjjTcu73/o0CEWL17MlClT6Ny5c/J+oCIiIi7WWiIjv+Xixd40bHg3c+bMIVeuXE7HSrF0XUw5ISwsjNDQUABCQ0MJCwtj4cKFvPjiiwCUK1eOcuXKXd4/PDycihUrUqFCBTZt2sTmzZsvb3v66acBKFu2LFWrViUgIIA8efKQJUsWoqOjqVy5Mr/++itffPEFGzZsICAg4Ia5nn/+ecaOHUtYWBjNmjW7vP7kyZNER0dTq1YtAFq1asXChQtvu/3Lli2jefPmALRs2ZLFixdf3ta4cWN8fHwoVaoUR1LaBSgiIl7j4kXo2zeedu1e57vvPuLllzczefIYsmXL5nQ0t0jXV/PNn///z/38rlxOSu7cV+6fnLFUJ06cYO7cuWzYsAFjDHFxcRhjqFChwnX33717Nz179mTlypUEBQXRunVrYmJiLm+/1Hvk4+NzRU+Sj48PsbGx1KxZk4ULFzJ16lRat27NO++8Q0BAwOWepcSDzvPly4efnx+zZ8/mp59+StY4p0yZMhEfHw9wRa7blbgNTt3TUUREPM/48efp1CkLsJ+PP/6Yr776CmOM07HcRj1TiYwbN46WLVuyd+9e9uzZw/79+ylSpAiVKlVi1KhRAGzcuPHyqbxTp06RPXt2cuTIwZEjR5g+ffotfd7evXvJmzcvr7zyCm3btmX16tU0adKEtWvXsnbtWkJCrryfYteuXfnuu++umA02R44cBAUFXR5v9dtvv13upSpcuPDlU5I3GuMUEBDA6dOnLy9Xq1bt8lixkSNHeuyVFSIikj5ERUXRr98jwIP06dOQ7t27Z6hCCpLRM2WMGQo8CRy11l4z1N4Y0wL4EDDAaeB1a+06dwe9Wt681x9snjfv7b9nWFgYH3744RXrnn32WdasWcO5c+coWbIkJUuWpFKlSgCUL1+eChUqUKJECQoUKMDDDz98S583f/58vv/+e/z8/LjjjjsYMWLETfevVq3addcPHz6c1157jbNnz3Lffffx66+/AvDee+/x/PPPM2jQIJ544onrvvapp56iadOmTJo0iZ9//pmff/6Zl19+me+//548efJcfi8REZFbsXcvhIae58SJ5uzZs4IxY37j+eefdzpWqjBJna4xxtQEzgAjblBMVQO2WGujjDGPAV9Ya6sm9cEhISH26rmatmzZQsmSJW8lv6RDOo4iIjJ16i4aN/Yhc+ZX+PPPj6lbt67TkVLEGLPKWhtyvW1J9kxZaxcaYwrfZHviwTvLgfy3GlBEREQyhn374ODBZbz00pPkzJmVGTOm3HDscUbh7gHobYBbGzgkIiIiGUJEBNSoEQf8QqFCOZk5cyb33Xef07FSnduKKWNMHRKKqeo32acd0A644QRd1toMNzDNm+gqPxER77V+/TBiY49SvvxBZs5cwl133eV0pDThlqv5jDHlgMFAI2vt8RvtZ60dZK0NsdaG5MmT55rtWbNm5fjx4/oH2UNZazl+/DhZs2Z1OoqIiKQRa6F/f0uXLt/Trt3L1K8/h0WLJnlNIQVu6JkyxhQE/gBaWmu3p+S98ufPz4EDBzh27FhKY4lDsmbNSv78GjYnIuIt1qyJp2NHi7X7adGiBUOHDiVz5sxOx0pTyZkaIQyoDeQ2xhwAugB+ANbaAcDnQC6gv+v0XOyNRrsnxc/PjyJFitzOS0VERCQNWQsXLpznu+9ewtqdvPNOHb7/fgQ+Pt43hWVyruZrlsT2tkBbtyUSERGRdO3ECXj22YucPPkea9aE07NnT959912nYzkmXd9ORkRERNKfQ4eOsHz5SS5cOMRvv/12+f613krFlIiIiCTLvn1w7twOnnqqAT4+x5k2LZwGDRo4HctxKqZEREQkSfv3Q5kyscTFjSdbttPMn/8/Kleu7HSsdEHFlIiIiCRp06aZxMQsJW/eucydu5SiRYs6HSndUDElIiIiNzRiBERHT+Ddd5+ndOnSzJgxg3z58jkdK11RMSUiIiLXdfQovPrqeWJiDlCnTg0mTJhAjhw5nI6V7njfZBAiIiJyU9ZCfHw8PXq8R0xMJZ59dinTp09XIXUD6pkSERGRy2JioHnzeI4cGcjSpT/QsWNHevfuja+vr9PR0i0VUyIiInLZv/+eYcGCHZw4sZnu3bvz0Ucf4brDidyAiikRERHh77/h4sWjNG36BNHR6xgyZAD/+c9/nI7lEVRMiYiIeLmzZ6Fq1YtERy8jLm4Tkyb9wZNPPul0LI+hYkpERMTLbd++llOnhuPjs5ZZs+bw0EMPOR3Jo6iYEhER8VJTp8LevSvp3LkegYGBzJw5k5IlSzody+OomBIREfFCFy5A27anOXLkBKVKFWTGjBnkz5/f6VgeScWUiIiIF7E24c9Bg/py+PCPPPjgA0ybtoigoCBng3kwFVMiIiJeIj4e3nzTsmbNPJYs6USjRo0ICwvD39/f6WgeTTOgi4iIeIm4uFhmzVrEkiVreOWVdowbN06FlBuomBIREcngTp2C3bvP8swzTdi+vTZdupxm4MABZMqkE1TuoJ+iiIhIBmYtNGhwkY0bd3PmzDR++aU/r732mtOxMhQVUyIiIhnY/v37OHjwK86f/5vx48fyzDPPOB0pw1ExJSIikgGtXAnLlu2lR4+HOXPmDLNnT6JWrVpOx8qQVEyJiIhkQO3bR7F69Tny5vVh0aJFlC1b1ulIGZaKKRERkQzEWpg0aSIbNnSgSJH7mDNnEYUKFXI6Voamq/lEREQyiJ49oWrVnTzzTFOCgwuwfPkEFVJpQD1TIiIiGYC1llmz/sfKlSdo2PAJxo0bRfbs2Z2O5RXUMyUiIuLBLlyAHTviaN++PbNnP8pLL81g8uRxKqTSkHqmREREPFibNnGMGxdNTMzvdO7cma+//hpjjNOxvEqSPVPGmKHGmKPGmI032F7CGLPMGHPeGPOe+yOKiIjI9URHR7NlSxtiYt6md++v+Oabb1RIOSA5PVPDgL7AiBtsPwG8ATR2UyYRERG5iV27YPz4aH7/vSZbt24lLGwEoaGhTsfyWkkWU9bahcaYwjfZfhQ4aox5wo25RERE5AY+/zyK0aPjyZbtFNOnT6devXpOR/JqaTpmyhjTDmgHULBgwbT8aBEREY8XHw8rVixnxowmBAUVYNasP6hYsaLTsbxeml7NZ60dZK0NsdaG5MmTJy0/WkRExKONGQPlyp2gTp0nCQrKTkREmAqpdEJX84mIiHiAJUvmsmmToVy50syaFU7evHmdjiQuKqZERETSKWthyxbLn3/24OefO1Ov3iNMmDCFgIAAp6NJIkkWU8aYMKA2kNsYcwDoAvgBWGsHGGPyAZHAnUC8MeYtoJS19lSqpRYREfECX30VT9euF4mNHUizZs0YNmwYmTNndjqWXCU5V/M1S2L7YSC/2xKJiIgI58+fZ9WqN4iNzcxbbzXmhx964uOjG5ekRzrNJyIiko788w/8/HMMixc/ydy5c+jRowfvvfeeJuNMx1RMiYiIpCMDB56iW7fM+PgcY8SIEbRs2dLpSJIEFVMiIiLpQHw8/PXXToYMaUCWLNmZMKEHDRs2dDqWJIOKKREREYctWACvvHKW48eb4ONzigULwqhSpYrTsSSZVEyJiIg4bNOm5eza5cvdd9/B3LlLKFasmNOR5BaomBIREXHIhg2wcWMYb73VitKlSzJjxnTuuecep2PJLVIxJSIi4oDRo6F5c4u1A6hd+2EmTpxIjhw5nI4lt0HFlIiISBqLj48nIuJTrL3AM8/kY+TI4WTNmtXpWHKbVEyJiIikkbNnoXv3OPbsaceoUUNp3749ffqMwtfX1+lokgIqpkRERNLIzJnn+OabzFh7kG7duvHJJ59oMs4MQMWUiIhIKouLg6iof/j22yeA4/z3v51p27at07HETVRMiYiIpKLNm6FRowtcuNCWo0fXM3HiGJ5++mmnY4kbqZgSERFJRQcObGbv3rNkzfoP//vf/3j44YedjiRupmJKREQkFaxdC1FR83nuuUbcddedzJw5g9KlSzsdS1KBiikRERE3i4iAhx6y+PiEUbx4fmbMmEGBAgWcjiWpRMWUiIiIm0VG9sfa7VSuvIupUxeRM2dOpyNJKlIxJSIi4gZxcfDll5YzZ76jV6+PePrppxk9+k/8/f2djiapTMWUiIiIG2zcGMvXX8cRF7eXtm3b8ssvv5Apk/6Z9QY+TgcQERHxZHFxcPbsWT7//Fni4orz6ae5GTRokAopL6IjLSIicpsOHYJHH40lNvZLtm37k379+tK+fXunY0kaUzElIiJym06ePMDu3fs4f34t4eHhNG3a1OlI4gAVUyIiIrdo/XqIj9/MU081wMfnFLNnT6J27dpOxxKHqJgSERG5BQcPQtWqcVg7m6CgWBYtWkj58uWdjiUOUjElIiJyC1atmkxs7BQKFtzAnDnLKFy4sNORxGEqpkRERJJgLfzwA/z77zi6dn2BkJAQpkyZQp48eZyOJumAiikREZEk/POP5csvz3DmzB4aNnyUcePGkT17dqdjSTqR5DxTxpihxpijxpiNN9hujDF9jDE7jTHrjTEV3R9TREQk7cXGQlxcHF980ZEzZ0rSsuVGJk+erEJKrpCcSTuHAQ1vsv0xoKjr0Q74JeWxUiZfPjDm2ke+fE4nExERT3HmDNSrF0fFir/Tv39/PvigBcOH/4qfn5/T0SSdSfI0n7V2oTGm8E12aQSMsNZaYLkxJtAYc7e19pCbMt6yI0dubb2IiMjVLlw4ydatKzh6dBa9evXirbfecjqSpFPuGDN1L7A/0fIB1zrHiikREZHbtWMHXLhwiObNGxIVtYVRo4bTrFkzp2NJOpamA9CNMe1IOBVIwYIF0/KjRUREkhQTAzVqXOTUqfX4+v7F1KlTqV+/vtOxJJ1zRzF1ECiQaDm/a901rLWDgEEAISEh1g2fLSIi4jbr16/g7NmfyJp1J7Nnz6dSpUpORxIPkJwB6EmZDLzkuqrvQeCkk+OlktKv3x9ORxARkXRmxAjo3n0lderUIU+e5axYMVKFlCRbkj1TxpgwoDaQ2xhzAOgC+AFYawcA04DHgZ3AWeDl1AqbXHnzXn+wua/vGTp2fJZdu96mR48eZMqkabZERLzdxYvw2Wcn2L//b4KDizN9+nTy5s3rdCzxIMm5mu+mo+5cV/F1cFsiNzh8+PrrL17MwrvvdqJXr1mMGPEcixYVpWTJ3GkbTkRE0oX4eIiPt/Tu3ZN9+3pQp04IEyfO584773Q6mngYd5zm8xh+fn706dOHjh1/5vjxfDzyyGOsXr3a6VgiIpLG4uOhRQtLhQoRfPDBB7zwQj2mT5+oQkpui1cVU5f8/HMdli49ga/vEapVe5hPP53pdCQREUlDsbEX2LBhIhs3TuCNN95k1KhRZMmSxelY4qG8spgCeOihSkRGRnLffZ3p3r0BTZoM4OLFi07HEhGRVHTkCKxbd4Ynn3ySTZue4dtvc9K7dy98fLz2n0NxA6/+7bnrrruIjPyYevXCmDixPY8++ihHjx51OpaIiKQCa+GJJy7y4IOHmTNnHsOGDePDDz/EGON0NPFwXl1MAWTL5sf//teMESOGs2zZdgoW/IvfftvsdCwREXGzv/7axeHDL2BtG/78cxKtWrVyOpJkEJobwKVly5b4+1eiWbMstGnTCmvb8dJLLzkdS0REUmj2bFiwYC///W814uLiWLBgKlWrVnU6lmQgXt8zlVjTpqXYt+9OqlfPTKtWrWjSpC9nz2oclYiIJ+va9SjffBNFlix3sGTJEhVS4nYqpq5y9915mDVrFq1bd2fixNcpVeo3jlxvBlAREUnXLlyA0aNHs3x5UUqUaM/y5YsoXry407EkA1IxdR2ZMmXi118/5q23FnP48AeEhISwYsVKp2OJiEgyffwxlC27j2bNWlGtWjBLlkzjnnvucTqWZFAqpm6iV69aLFs2Gx+fLDz44Anatl3kdCQREUmCtZb168PZvv1PmjR5mpkzZxIYGOh0LMnAVEwloUKFCsybF0FgYB6GDBlEx44duXDhgtOxREQEyJcPjLny4eNjmDq1Jq+9tpGxY0eTNWtWp2NKBqdiKhnuuy8Xhw+X491389KvXz8qV36LDRuOOR1LRMTr3XhIaz769++Pr69vWsYRL6ViKpkyZ85Ez549GT58NOvXf0LlyuuJiIhwOpaIiNyAJuOUtKJi6ha99NILjBlzhty5v6BmzZr8979DnI4kIuKV4uKcTiCSQMXUbXj++eKsWzeRmjVr0q5dDKVLL+T8eY2jEhFJC2fPJvy5bt1qZ4OIuKiYuk25cuVi2rTpVKlSis2bl1GvXl0OHz7sdCwRkQxt3jwoUCCe5577kpCQEKfjiAAqplLEzy8TERF1CAsrzJo1ayhf/nkGD97odCwRkQwnJgZiY2OJiBjE6dOT+OOPMN58803uvPPsdffPkeNcGicUb6Ziyg1CQ19g2bJlnD79Ba+8koP+/TWOSkTEXdq2herVTxAcXIGPPnqVWrX6s2HDH/Tq1YuTJ7Px++8jKVSoMMb4UKhQYX7/fSTR0f5OxxYvohsdu0m5cuXYsCGKF1/8lA4d+rNu3Up+/PEnsmfP4nQ0ERGPc/48ZMkCe/fuZc2aGaxevY9Chc4yYcIEGjVqdMWVei1atKBFixYOphVvp54pN7r//iAWL+5D586dGTTIki/fdjZvPuR0LBERj7JtGxQtGk/z5iMpUaIEW7a8Tbdu/mzZspHGjRtrygNJd1RMuZmvry/ffPMN77zTinPntlO3bghLly51OpaISLp34ULCrWBWrx7PP/9MJyysL40aNWLr1q18+umn+Pvr1J2kTyqmUskPP1RjzZpi3HGHP7VqNeLll+dirdOpRETSpx49oEyZGOrUeZTmzZvywAMfMX/+t4wePZqCBQs6HU/kpjRmKhWVLVuWlStXUqXKnwwbVoOzZz9lxIjPyJJF46hERC5NunnqVBSLFo1h5044dmwb/fr1o127dmTKpH+ixDOoZyqVBQUFsWlTC156aTDh4d2pVasWf/110OlYIiKOioqCqlUtL764hKJFizJtWgdee209O3euoX379iqkxKOomEoDmTP7Mnz464wdO5Z163JQtCgMG7bG6VgiImkuNjbhz02blrB79xRGj/6BUqVKsWrVKvr370+uXLmcDShyG1RMpaGmTZsyenQ/smbdTtu2j/DLL79gNZBKRLxEeDjcf38szz3Xjho1quPv/zphYc+zYMECgoODnY4nctuSVUwZYxoaY7YZY3YaYzpfZ3shY8wcY8x6Y8x8Y0x+90fNGBo1eoCDByvQoMGDtG/fkYceGk10dIzTsUREUk1sLJw/f54lS4Zy8OB0Jk2azSeffMK2bdsIDQ3VVAfi8ZI8KW2M8QX6AfWBA8BKY8xka+3mRLv1BEZYa4cbY+oC3wAtUyNwRhAYGMjkyZN56aXfGDWqGZUrd2bevI7kz68aVEQyjrg4eOYZi4/PX2zY0IBdu3bRqFEjfvjhf9x///1OxxNxm+T0TFUBdlpr/7LWXgBGA42u2qcUMNf1fN51tstVfH19GTmyNd9/P4fDh/tRqVIlZs1a4nQsEZEUu3SV3s6d21i/fiITJ/bFz8+PmTNnMnHiRBVSkuEkp5i6F9ifaPmAa11i64BnXM+bAAHGGI0iTIb33qtHREQE2bKVpEGDIrz88iyNoxIRj7V0KTzwQDyvvNKDsmXLcuJEa378sSDr16/n0UcfdTqeSKpw1wD094Baxpg1QC3gIBB39U7GmHbGmEhjTOSxY8fc9NGer1SpUsyfP5F77tnOsGEdaNOmDTExGkclIp4jPh7i4+OJjAzn77+XMXjw77Rs2ZLt27fz9ttv4+fn53REkVSTnIk8DgIFEi3nd627zFr7N66eKWPMHcCz1troq9/IWjsIGAQQEsOCknkAAB1iSURBVBKi7pdEChUKZP/+mnzxRTO6devGnDklCA9/kapV73E6mojITb33HqxZc4J//32ciIgIqlatSp8+g6lSpYrT0UTSRHJ6plYCRY0xRYwxmYFQYHLiHYwxuY0xl97rI2Coe2N6Bx8fH7p27cqQIdPZt+816tYdx8KFC52OJSJyjfj4hD+PHDnCwoXhzJ37G7t372PYsGEsXbpUhZR4lSSLKWttLNARmAlsAcKttZuMMV2NMU+7dqsNbDPGbAfyAt1TKa9X+M9/GjJ9+lHy5x9EvXr16NFjIPHx6sgTkfRh1y4IDrZ06DCWYsWKsWZNC9577wA7dmylVatW+PhoCkPxLsapwc4hISE2MjLSkc/2FCdPnqRZszZMn/45JUocZvXqGrpruog4xlowBqZMmcMLLwRw9uwnNGjgS+/evSlRooTT8URSlTFmlbU25Hrb9N+HdCxHjhxMnhxO7don2Lr1B2rUqMG+ffucjiUiXuiXX6Bq1XM0bvwsTz31CPnyNWPy5DeYPn26Cinxeiqm0rlMmXyYN682kyZ1YMeOHZQp043evdc6HUtEvIC1CY9///2X2bPDWbVqOrNmLeXrr79m06ZNPPXUU5q9XITkXc0n6cDTTz/NsmUrqFgxjrffPgD8xJtvvqG/yEQkVURHQ2iopUiRFUyZ0pQDBw7QvHlzvvtupe7WIHIV9Ux5kFKlirNzZ34ee+w33n77LZo3b8eJE+ecjiUiGcilYbS7d69jxYoVDBgwgDx58rBo0SJGjhypQkrkOlRMeZj8+e9kypThdO3aldGj61OgwC527drrdCwRyQCmToXg4FjatXuHkJCKGPMEAwc+xMqVK6levbrT8UTSLRVTHsjHx4fPPvuMzz8vTHz8UKpWrcTcuXOTfqGIyHVYC3FxcSxY8AebN69iyJBJdOjQgZ07d9CuXTt8fX2djiiSrmlqBA+3fft2GjduzLZtgTz11Nf88UctfHw0jkpEkhYXB+3agbV7WL26MevWraN27Tr06fMTZcuWdTqeSLqiqREysGLFihEREUGhQp8zaVJBQkNf4ezZs07HEpF07NL/of/+ez9z5y7g119HEBUVxdixY5k7d44KKZFbpGIqAwgICGD79kf58MM/GTduKNWqVWfFiv1OxxKRdGj1aggJief99/tQokQJDh1qSJcu8WzZsoWmTZvqCmGR26BiKoPIlMmHb799kylTprB1ayMefDAbo0YtcjqWiKQj1loiI2exYcMOevb8nccee4ytW7fwxRdfkC1bNqfjiXgszTOVwTz++OP8+WcpWrQI48UX3+TQoR688847+t+miBfr3h3WrYsiOvoFZs+eTalSpenT5yfq1avndDSRDEE9UxlQ/fqF2bWrFc8804T33vuR4sVnc+yYxlGJeKOTJ08yefJsxo2bQUREJD/99BNr165RISXiRiqmMqiAgADGjh3Lc88NYMeOqtSo0Zzdu3c7HUtE0sC+ffDII5bPPptMsWLFWLGiAW3bzmPnzm288cYb+Pn5OR1RJENRMZWBGWMID3+KsLCVHDmygJCQEEaMWOx0LBFJZTt2rGTJkr/46quh3H///URGrmTQoEHkyZPH6WgiGZKKKS8QGvoIK1euJCDgeVq1epA2bcJxan4xEUkdI0fC44/H0KpVax55pAqBgTX5/ffnWLJkCZUqVXI6nkiGpgHoXuKBBx4gIuJ76tSZyNChLTlzZjxDhw4le/bsTkcTkRS6cOECEybMZdasO/DxmU7nzp35+OOPCQgIcDqaiFdQMeVF8ua9g02bnqVHj1107vwls2bNYfz48tStW8jpaCJyi06dgg4doGDBNYwbF8r27Tt44okn6NVrMUWLFnU6nohXUTHlZYwxfPjhh9x5Zw06dCjO0093Yvz4VjRo0MDpaCJyC/7+eydTpsQRHT2MokUtU6dO4fHHH3c6lohX0pgpL/X669VYvfok9923kccff5x33x1MfLzGUYmkZ/PnQ/36sbz//meUL1+aixer0KNHfjZu3KhCSsRBKqa8WHDwfSxbtoyGDV/nxx+bUbZsOGfOnHE6lohch7WWmTPnMH/+Hnr2HE1oaCg7dmzl/fffJ3PmzE7HE/FqxqmrukJCQmxkZKQjny1Xio+3PPPMHCZPfonSpXMxYcIEHnjgAadjiXi9uDj45BOIjd3P8uXNWLJkCRUrVqFv39489NBDTscT8SrGmFXW2pDrbVPPlODjY5g48RFmzhzOwYN/U7r0fLp1W+l0LBGvd+LEMUaPXssPP0xi+/btDB48mJUrl6mQEklnVEzJZfXr12fevNX4+FTn888n8fXXX2s+KpE0tnkzPPZYPF9/PZhixYqxf39V3n77L7Zv306bNm3w8dFf2yLpjb6VcoXy5Quxf38BQkP/4pNPPqFBgzf4++/TTscS8RorViznf/87yiefDCckJIQNG9bw448/EhgY6HQ0EbkBTY0g18idOzujRo2kYsUqfPDBUxQrtpY1a/Jp7hqRVPLLL7BhQzRHj7Zl/PjxFCr0AL17f0+jRo0wxjgdT0SSoJ4puS5jDO+//xbff3+STJm+pHLlykydOtXpWCIZztmzZxk2bDkDB65g6tSZdOvWjS1b1tO4cWMVUiIeQsWU3NR771Vk7drBFClShCefXMAjj8wnLi7e6VgiHu3wYWjWzPLjj7MoWbIkK1bUpGnToWzfvplPP/0Uf39/pyOKyC1IVjFljGlojNlmjNlpjOl8ne0FjTHzjDFrjDHrjTGaPS4DKVy4MIsXL6FIkaeYM+cITZs25fRpjaMSuV07d25hwoR/ePfd4QQGBjJ//mzGjBlNgQIFnI4mIrchyWLKGOML9AMeA0oBzYwxpa7a7VMg3FpbAQgF+rs7qDgre/Zs7NxZnR49DvPnn5OpVOlxZsz4y+lYIh5j6lRo3fo8nTq9Qe3aZfH3L0u/fg+zatUqatWq5XQ8EUmB5AxArwLstNb+BWCMGQ00AjYn2scCd7qe5wD+dmdISR98fAzvv/8mlSqV5bHHzvL449kYO3Yazz6rjkiRm4mLi2P48Aj++COQ+PhRvP76q3Tt2pVcuXI5HU1E3CA5p/nuBfYnWj7gWpfYF8CLxpgDwDSg0/XeyBjTzhgTaYyJPHbs2G3ElfSgbt26zJ0bTJEiX9O06RN07dqV+HiNoxJJ7N9/4Z134McfN1G5cmXGjq1JtWodWbNmDv369VMhJZKBuGsAejNgmLU2P/A48Jsx5pr3ttYOstaGWGtD8uTJ46aPFic8/HB+Nm78jhdffJEuXSLInz+CgwdPOR1LJN04evQgQ4Yc5N13wzh27BijR49kwYI5lC9f3uloIuJmySmmDgKJR0Xmd61LrA0QDmCtXQZkBXK7I6CkX/7+/owYMYLnn3+XQ4d8qVu3Flu3bnU6lohjVq2C5s3j+Oqr7yhbtjgxMaX59FPD1q1beeGFFzTVgUgGlZxiaiVQ1BhTxBiTmYQB5pOv2mcfUA/AGFOShGJK5/G8gDGGMWPq8r//nePEiQNUrvwwX321xOlYImnOWssff6wgPPwfPvtsGPXr12fLltV069aN7NmzOx1PRFJRkgPQrbWxxpiOwEzAFxhqrd1kjOkKRFprJwPvAv81xrxNwmD01lY3dfMq9erVYtWqVVSrNonPPnuIQ4f68vPP7XUfMcnQ4uLgp58gJuYQixb9hxkzZlCsWDA///wTjz76qNPxRCSNGKdqnpCQEBsZGenIZ0vqOXHiHM88M5gFC97gqaeeYsSI3wgMzOF0LJFUcfLkKcqWPcqBAxEEBLTniy++oGPHjvj5+TkdTUTczBizylobcr1t6jYQt8qZ05958zrSp08fpk79i3z59jJ9+i6nY4m4zZ490KqVpX//URQvXoz9+4N5+eW5bN++nbfffluFlIgXUjElbmeMoVOnTvz882/ExmbiueeeZOLEiU7HEnGLxYs3MXLkv3ToMITChQuzYsU8hgwZQt68eZ2OJiIOUTElqaZ9+wrs2hVAqVIBNGnShOefDyM2VvNRiecJD4fPPz9NmzZtaNmyDDlzlmf48FYsXbqUypUrOx1PRByWnBnQRW5boUIFWLhwIY0b92Hs2GZs3/4j8+f/h8DAQKejiSTLxYsX+eGHbaxadQ5jRvH+++/z6aefcueddyb9YhHxCuqZklSXNWtWpk17n1dfncbGjZ2pUqUKGzduTvqFIg6JioI33oDhwxdTvnx5Vqx4kEce+YKNG9fSo0cPFVIicgUVU5ImfHwMAwY8zrx5c4iKgvLlL/LRR0udjiVyXVu37mXAgDO0bj2MCxcu8Oefo5k+fQrFixd3OpqIpEMqpiRN1ahRgxkz5pM9ezzffvsen376KXFxcU7HEmHJEnj//Qt89tln1KlTHD+/4nzzzQNs2rSJJ598UrOXi8gNacyUpLlKle7hyJGcdOxYku7duzN9uh9//PEWhQppPipxhrWWPn02Mn58TuLi+tGiRVO+++477r336nu6i4hcS8WUOMLfPyuDBw+mePEafPjhc5QvP57FiytSpkwZp6OJlzh/Hn78Ee6+eye//tqGhQsjKF++In37TqZ69epOxxMRD6LTfOIYYwwffNCagQN3kjlzVx588EHGjBnvdCzxEkePHufrr0/w8svj2bRpEwMH9mHVqkUqpETklqmYEse1a1eeNWsWUKZMMKGhAdSoMUfjqCRVbN8Ob70VT9++/Slfvihnz5agU6eD7Nixg3bt2uHr6+t0RBHxQCqmJF249957mT17DiVKZGHx4lE8+eSTREVFOR1LMpjfftvCzz//S6dO/QgODmbdurn06dOHoKAgp6OJiAdTMSXpRkBAFrZsqcXAgVWZM2cOpUu3Z/z47U7HEg9mLfz2Gwwc+A+hoaF89VUp7r67BuPGdWXOnDkaoycibmGstY58cEhIiI2MjHTksyX9W7x4KXXqBBIff4ZRo/bwwgvPOx1J0rF8+eDIkWvX33WXJXPmv/n779Vkzvw8nTt35v333ydbtmxpH1JEPJoxZpW1NuR629QzJelS9erViIjISfnyPQgNfYEPPviYCxc0jkqu73qFFMDRo4YDByrQpMlvbN26lS5duqiQEhG3UzEl6VbFivlYtmwkr776Kt9/n4d77lnD4cMnnI4lDrAWYmP///m8eQmDyQEuXrz5a+fMGc24ceEUKlQodUOKiNfSPFOSrmXJkoUBAwZw9uxCRo5cSLVqzzNx4kTKlSvndDRJgX37Ev4sWDDhz2HDEk7VNWyYsNy8OVStamnT5l+io6MpXTofjz9+mObN1xAVFc1//tOchx5aTpUqfxAVFQ0MueFn1a1bN1XbIiKiMVPiMZYvX84zzzxDVFQuXn11AL17P+x0JK8VGwv//gs5XJPWL1mSMAnmpbrl++/BGEuHDjFER0fTrNkdZMt2jo4dI4mOjuaddx4jT57DNGw4lOjoaEaP7kr27FsoUOBDoqOj2bu3P7Gx87H2a9cnfgSsAWa4lqsB+8iePYrAwEAOHjxww6wO/RUnIhnMzcZMqZgSj3Lo0CGCgyM5erQaHTv2plevLmTKpA7WlNq3L2HcUeXKCcuTJsHu3fD66+eJjo6ma1dfdu6Ed95ZRXR0NN26PUR0tB8vvvgT0dHRTJr0OufOZaZEiZeJiopiz55eXLx4Dmubuj7hc+Ak8JNr+XEgiqxZ1xAYGEhAQCFy5sxMzpx3EBgYeMUjKCjomnWXHn5+fgDc7LZ5KqZExB1UTEmGcubMBdq0+Y7w8M955JFHGDp0NAUK5HI6lqPOnYNjx6BAgYTCYv16WLECWrW6yMmTJ/n991hmzPDjvffWEBUVxa+/PsDSpQ/QocO3REdHM2vWM+zbV5FKlZ4gOjqaPXs+JSamBta6zsPxMXA/0Ma13Bi4Az+/MQQFBZEtWzly5MhO3rznkiyCLq3PkSMHWbNmdUv7AwPPcfKk/zXrc+Q4R3T0tetFRG6ViinJkIYMGcKrry7HmM+ZPPkUjz1W2ulIKRIXl1AI+fgk9BKtWQMPPxxHbOwpFiw4y9ixmQgN3crFi8eZOTOIiRPL0arVAM6dO8ySJQ+zZk0oDz1Un9OnD7Nv34ucOvUhkBU4D7xJQiEUDMQDjYBq+Pp+TGBgIP7+VfD3L0jBgjsJCgoiICAXOXMGEBR08x4if39/zM26hdLQyJEj+eSTT9i3bx8FCxake/futGjRwulYIpJBqJiSDGvo0I107LgVa19myJCBNG/e3OlIAMTEwJ49cPfd8cBptmw5zfjxhipV9pE581FWr/YlPLwM9er9SaZMu1i37j7mz+9ISEgbLlxYzcGDNTl+/GegBLANeBHoB5QH9gCPklActScwMI7s2R/C17cqhQsvI1eubPj730v27Lm5914fcua8/imyoKAgsmfPnm6KIRGR9EzFlGRohw8f5rnnnmPx4nFA3mu2580Lhw/f2ntam3DqDCzx8f9y6NBJpk2L5e67j+Lvf5g9e2IYPbokpUtHEBCwmT17/Jkx4w0eeOBHfH1ncfhwIQ4dmoAxT2LtVBIGTC8hoQiaDVQGfgZe5c47d5M9e0Xi4kIpUmQWefPG4u9fEGvv54EH/iVPnuzkyBFIUNC1vUMBAQH4+GiGExGR1KZiSjK8CxcukCVL5htuP3jwEEuXxmDMCfz9j3D8+EnGji1C3rw7CQpaz9Gj/zJ1agfy5p2Mv/94/vnHhz17IvDxeZP4+D5APuAQ8BowELgXiAQ6kT37dAICihIT04VChaZRoMAh/P3v4eTJhyhZ8igFCvgSEJCLrFlzcc892a/oKcqRI4durisi4gFuVkzpMijJEDJnvnEhBXDvvfcA54BwoLNr7UVgFlmz/kxgYC7Onn2JCxfOUrBgHu6/Pxf58k2gVKlilCjRg4CAIE6fnsUDD7xA/vxtXD1EmciRY9TlK8oSNE6V9omISPqlYkq8wi+//MKePREULVqNMmWWuwZPH+euuz4ga9bPE+1Z1bGMIiLimZJVTBljGpIwQYwvMNha++1V23sBdVyL2YC7rLWB7gwqkhKvvfaa0xFERCSDSrKYMsb4knAZUX3gALDSGDPZWrv50j7W2rcT7d8JqJAKWUVERETSneRcBlQF2Gmt/ctaewEYTcIkNTfSDAhzRziRW5H32gv5brpeRETEHZJzmu9eYH+i5QPcYGCJMaYQUASYm/JoIrfmVqc/EBERcQd3T1ATCoyz1sZdb6Mxpp0xJtIYE3ns2DE3f7SIiIhI2ktOMXUQKJBoOb9r3fWEcpNTfNbaQdbaEGttSJ48eZKfUkRERCSdSk4xtRIoaowpYozJTELBNPnqnYwxJYAgYJl7I4qIiIikX0kWU9baWKAjMBPYAoRbazcZY7oaY55OtGsoMNo6NaW6iIiIiAOSNc+UtXYaMO2qdZ9ftfyF+2KJiIiIeAbH7s1njDkG7E2Dj8oN/JMGn5Meqe3ey5vb781tB+9uv9ruvdKi/YWstdcd8O1YMZVWjDGRN7oxYUantntn28G72+/NbQfvbr/a7p1tB+fb7+6pEURERES8ioopERERkRTwhmJqkNMBHKS2ey9vbr83tx28u/1qu/dytP0ZfsyUiIiISGryhp4pERERkVTjscWUMaahMWabMWanMabzdbZnMcaMcW2PMMYUTrTtI9f6bcaYBmmZ2x2S0fZ3jDGbjTHrjTFzXDegvrQtzhiz1vW4ZiZ7T5CM9rc2xhxL1M62iba1MsbscD1apW3ylEtG23slavd2Y0x0om0efeyNMUONMUeNMRtvsN0YY/q4fjbrjTEVE23z6OMOyWp/C1e7Nxhjlhpjyifatse1fq0xJjLtUrtHMtpe2xhzMtHv9+eJtt30O5PeJaPt7ydq90bX9zyna5tHH3cAY0wBY8w8179pm4wxb15nH+e/+9Zaj3sAvsAu4D4gM7AOKHXVPu2BAa7nocAY1/NSrv2zAEVc7+PrdJvc3PY6QDbX89cvtd21fMbpNqRB+1sDfa/z2pzAX64/g1zPg5xukzvbftX+nYChGejY1wQqAhtvsP1xYDpggAeBiIxw3G+h/dUutQt47FL7Xct7gNxOtyEV214bmHKd9bf0nUmPj6TaftW+TwFzM8pxd7XhbqCi63kAsP06f+c7/t331J6pKsBOa+1f1toLwGig0VX7NAKGu56PA+oZY4xr/Whr7Xlr7W5gp+v9PEWSbbfWzrPWnnUtLifh5tQZRXKO/Y00AGZba09Ya6OA2UDDVMqZGm617c24yY3HPY21diFw4ia7NAJG2ATLgUBjzN14/nEHkm6/tXapq32Qwb73yTj2N5KSvy/ShVtse4b6zgNYaw9Za1e7np8m4bZ29161m+PffU8tpu4F9idaPsC1P9zL+9iE+wueBHIl87Xp2a3mb0NCxX5JVmNMpDFmuTGmcWoETGXJbf+zru7eccaYArf42vQq2fldp3aLAHMTrfb0Y5+UG/18PP24346rv/cWmGWMWWWMaedQptT2kDFmnTFmujGmtGud1xx7Y0w2EgqF8YlWZ6jjbhKG61QAIq7a5Ph3P1n35hPPZIx5EQgBaiVaXchae9AYcx8w1xizwVq7y5mEqeZPIMxae94Y8yoJPZR1Hc6U1kKBcdbauETrvOHYez1jTB0SiqnqiVZXdx37u4DZxpitrh6PjGI1Cb/fZ4wxjwMTgaIOZ0prTwFLrLWJe7EyzHE3xtxBQqH4lrX2lNN5ruapPVMHgQKJlvO71l13H2NMJiAHcDyZr03PkpXfGPMI8AnwtLX2/KX11tqDrj//AuaTUOV7kiTbb609nqjNg4FKyX1tOncr+UO5qrs/Axz7pNzo5+Ppxz3ZjDHlSPidb2StPX5pfaJjfxSYgGcNbUiStfaUtfaM6/k0wM8YkxsvOvbc/Dvv0cfdGONHQiE10lr7x3V2cf6778SAspQ+SOhR+4uE0xiXBhWWvmqfDlw5AD3c9bw0Vw5A/wvPGoCenLZXIGHQZdGr1gcBWVzPcwM78LzBmMlp/92JnjcBlrue5wR2u34OQa7nOZ1ukzvb7tqvBAkDT01GOvau7IW58SDkJ7hyEOqKjHDcb6H9BUkYA1rtqvXZgYBEz5cCDZ1ui5vbnu/S7zsJBcM+1+9Bsr4z6f1xs7a7tucgYVxV9gx43A0wAuh9k30c/+575Gk+a22sMaYjMJOEqzWGWms3GWO6ApHW2snAEOA3Y8xOEn7JQl2v3WSMCQc2A7FAB3vlqZB0LZlt/x64AxibMOaefdbap4GSwEBjTDwJvZLfWms3O9KQ25TM9r9hjHmahON7goSr+7DWnjDGdANWut6uq72ySzxdS2bbIeF3fbR1/W3i4vHH3hgTRsJVW7mNMQeALoAfgLV2ADCNhKt6dgJngZdd2zz6uF+SjPZ/TsK40P6u732sTbjxa15ggmtdJmCUtXZGmjcgBZLR9qbA68aYWOAcEOr6/b/ud8aBJty2ZLQdEv7TOMta+2+il3r8cXd5GGgJbDDGrHWt+5iE/zykm+++ZkAXERERSQFPHTMlIiIiki6omBIRERFJARVTIiIiIimgYkpEREQkBVRMiYiIiKSAiikRSfeMMYHGmPau5/cYY8Y5nUlE5BJNjSAi6Z7rnlxTrLVlHI4iInINj5y0U0S8zrfA/a5J+3YAJa21ZYwxrYHGJMzwXBToScJM1y2B88Djron77gf6AXlImNTvFWvt1rRvhohkRDrNJyKeoDOwy1obDLx/1bYywDNAZaA7cNZaWwFYBrzk2mcQ0MlaWwl4D+ifJqlFxCuoZ0pEPN08a+1p4LQx5iTwp2v9BqCc627z1fj/2ytBwr05RUTcQsWUiHi684mexydajifh7zgfINrVqyUi4nY6zScinuA0EHA7L7TWngJ2G2OeAzAJyrsznIh4NxVTIpLuWWuPA0uMMRuB72/jLVoAbYwx64BNQCN35hMR76apEURERERSQD1TIiIiIimgYkpEREQkBVRMiYiIiKSAiikRERGRFFAxJSIiIpICKqZEREREUkDFlIiIiEgKqJgSERERSYH/A/F/va2B08XbAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "gmxzs4CQb5TF" }, "source": [ "## Local Error \n", "The Error for the 2 step Adams Moulton is:\n", "\\begin{equation}y_{n+1}=y_n+\\frac{h}{12}[5f(t_{n+1},w_{n+1})+8f(t_{n},w_{n})-f(t_{n-1},w_{n-1})] +\\frac{-h^4}{24}y^{(4)}(\\eta),\\end{equation}\n", "where $\\eta \\in [t_{n-1},t_{n+1}]$.\n", "\n", "Rearranging the equations gives \n", "\\begin{equation}\\frac{y_{n+1}-y_{n}}{h}=\\frac{1}{12}[5f(t_{n+1},w_{n+1})+8f(t_{n},w_{n})-f(t_{n-1},w_{n-1})] +\\frac{-h^3}{24}y^{(4)}(\\eta),\\end{equation}\n", "For our specific initial value problem the error is of the form:\n", "\\begin{equation}\\frac{-h^3}{24}y'''(\\eta)=\\frac{h^3}{24}2e^{-\\eta} \\leq\\frac{(0.3)^3}{24} 2\\leq 0.01042.\\end{equation}" ] }, { "cell_type": "code", "metadata": { "id": "4GouWupKb5TF", "colab": { "base_uri": "https://localhost:8080/", "height": 203 }, "outputId": "5c9e501c-4f0a-4f3d-b064-63e93bd1e7a0" }, "source": [ "\n", "\n", "d = {'time t_i': t, 'Adams Bashforth w_i': w,'Exact y(t_i)':y,'Error |y(t_i)-w_i|':np.abs(y-w),'LTE':round(2*0.5**3/24,5)}\n", "df = pd.DataFrame(data=d)\n", "df" ], "execution_count": 10, "outputs": [ { "output_type": "execute_result", "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
time t_iAdams Bashforth w_iExact y(t_i)Error |y(t_i)-w_i|LTE
00.01.0000001.0000000.0000000.01042
10.50.7130610.7130610.0000000.01042
21.00.7382410.7357590.0024820.01042
31.50.9491350.9462600.0028750.01042
42.01.2732551.2706710.0025850.01042
\n", "
" ], "text/plain": [ " time t_i Adams Bashforth w_i Exact y(t_i) Error |y(t_i)-w_i| LTE\n", "0 0.0 1.000000 1.000000 0.000000 0.01042\n", "1 0.5 0.713061 0.713061 0.000000 0.01042\n", "2 1.0 0.738241 0.735759 0.002482 0.01042\n", "3 1.5 0.949135 0.946260 0.002875 0.01042\n", "4 2.0 1.273255 1.270671 0.002585 0.01042" ] }, "metadata": {}, "execution_count": 10 } ] }, { "cell_type": "code", "metadata": { "id": "cp3BAQi5b5TH" }, "source": [ "" ], "execution_count": 10, "outputs": [] } ] }